HYDRO-UNIVERSITY COMPUTING CENTRE 
ALGOL PROCEDURES REFERENCE MANUAL 

This document contains details of Algol procedures which have 
been written for and tested on the Elliott 503. Tapes for these 
procedures are available in the Computing Centre. Users wishing to 
incorporate these procedures in programs submitted for punching 
should indicate their requirements on the Algol program sheet thus: 

TITLE J 

begin integer ; real ; 

LIBRARY LE04 



LIBRARY NC01 



comment then follows the program text; 



end ; 

Such instructions will be sufficient to ensure that the appropriate 
procedures are copied at the places indicated. It will be left to 
programmers to ensure that the order of copying of requested procedures 
is a valid order of declaration. This is important in those cases where 
one procedure uses another as indicated in the procedure descriptions. 

The manual also contains an index of HUCC Library procedures, 
together with an index of other known Algol procedures and their 
sources. In most cases copies of these are available in the Computing 
Centre. Procedure tapes , however, are available only for HUCC procedures. 



CLASSIFICATIONS 



C COMMERCIAL 

CA INTEREST 

D INFORMATION HANDLING-DATA PROCESSING 

DH SORTING 

DP PLOTTING 

DZ MISCELLANEOUS 

F FUNCTIONS-EVALUATION OF-ETC 

FB BESSEL 

FC COMPLEX 

FD POWERS AND EXPONENTIALS 

FE ELLIPTIC INTEGRALS 

FP POLYNOMI ALS-INC CHEBYSHEV ETC 

FS SPECIAL 

FT TRIGONOMETRICAL FUNCTIONS 

FZ MISCELLANEOUS 

G DIFFERENTIAL EQUATIONS 

GA ORDINARY-NOT LINEAR OR 1ST ORDER 

GL LINEAR 

GP PARTIAL 

L LINEAR ALGEBRA 

LA CHANGE FORM OF MATRIX 

LB BOOLEAN MATRICES 

LE LINEAR EQUATIONS AND INVERSION 

LF FORM SPECIAL MATRIX 

LG ARITHMETIC FUNCTIONS-ONE MATRIX 

LH ARITHMETIC FUNCTIONS-TWO MATRICES 

LL LATENT ROOTS 

LO DETERMINANTS 

LR READ OR INPUT 

LZ MISCELLANEOUS 

M MATHEMATICAL METHODS 

MC CURVE AND SURFACE FITTING 

MD DIFFERENTIATION-MAX AND MIN 

ME ERROR ANALYSIS 

MG GEOMETRY 

MH HARMONIC ANALYSIS 

MR ROOTS OF EQUATIONS 

MS INTEGRATION AND SUMMATION 

MT INTERPOLATION AND DIFFERENCES 

MZ MISCELLANEOUS 

N INTEGERS AND NUMBER THEORY 

NC PERMUTATIONS AND COMBINATIONS 

NP PRIME NUMBERS 

NR PARTITIONS 

NZ MISCELLANEOUS 

O OPERATIONAL RESEARCH 

OL LINEAR PROGRAMMING 

OP PERT-CRITICAL PATH ANALYSIS 

P PHYSICS 

PH HEAT 

PN NUCLEAR ENGINEERING 

PQ QUANTUM MECHANICS 

S STATISTICS 

SA ANALYSIS 

SB SMOOTHING 

SC CORRELATION 

SM MOMENTS 

SR REGRESSION 

SS STOCHASTIC PROCESSES 

SV ANALYSIS OF VARIANCE 

Z MISCELLANEOUS 

ZA COMPILER TECHNIQUES 

ZZ MISCELLANEOUS 



ABBREVIATIONS 



ACM COMMUNICATIONS OF ACM 

AIS ALGOL INFORMATION SHEET 

ALY ALGORYTMY 

BIT NORDISK TIDSKRIFT FOR INF BEHANDLG 

CJ COMPUTER JOURNAL 

JCM JOURNAL OF THE A C M 

MCA MATHEMATISCH CENTRUM AMSTERDAM 

NUM HUMERI SCHE MATHEMATIK 

ORN OAK RIDGE NATIONAL LABORATORY 

SUC STANFORD UNIVERSITY CALIFORNIA 

B001 SELECTED NUMERICAL METHS BY C GRAM 

NPL NATIONAL PHYSICAL LABORATORY 



INDEX of KUCC ALGOL MBMEY PROCEDURES 



DH 01 

DH 02 
DH 03 
DH 04 



DP 01 
DP 02 



FC 01 
FC 02 
FC 03 
FC 04 
FC 05 
FC 06 
FC 07 
FC 08 
FC 09 
FC 10 
FC 11 
FC 12 
FC 13 
FC 14 
FC 15 
FC 16 
FC 17 
FC 18 
FC 19 
FC 20 
FC 21 
FC 22 
FC 23 
FC 24 
FC 25 
FC 26 



FP 01 



FS 01 



GL 01 



LE 01 
LE 02 
LE 3 
LE 04 
LE 05 
LE 06 
LE 07 
LE 08 
LE 09 
LE 10 



SORT REAL NUMBERS INTO ASCENDING ORDER 

LOCATE ELEMENT IN LIST 

SORT REAL NUMBERS IN ASCENDING ORDER 

SORT ROWS OR COLUMNS OF MATRIX 



OUTPUT GRAPH OF VECTOR ELEMENTS 
OUTPUT PLOT OF ELEMENTS OF TWO VECTORS 



COMPLEX ASSIGNMENT 
MULTIPLE COMPLEX ASSIGNMENT 
COMPLEX MULTIPLICATION 
COMPLEX DIVISION 

ASSIMILATE REAL IN COMPLEX OPERATION 
ASSIMILATE IMAGINARY IN COMPLEX OPERATION 

COMPLEX CONJUGATE 
MODULUS OF COMPLEX NUMBER 
ARGUMENT OF COMPLEX NUMBER 
POLAR FORM 

MULTIPLICATION BY IMAGINARY OPERATOR 

COMPLEX SQUARE 

COMPLEX RECIPROCAL 

COMPLEX ROOT 

COMPLEX LOGARITHM 

COMPLEX EXPONENTIAL 

HYPERBOLIC FUNCTIONS 

COMPLEX SINE 

COMPLEX COSINE 

COMPLEX TANGENT 

COMPLEX INVERSE SINE 

COMPLEX INVERSE COSINE 

COMPLEX INVERSE TANGENT 

COMPLEX FQ9ER OF COMPLEX VARIABLE 

TEST FOR EVEN INTEGER 

LOGARITHM OF COMPLEX NUMBER 



SUM SERIES OF CHEBYSHEV POLYNOMIALS 



LOGARITHM CP FACTOR I AL ( n ) 



FOURTH ORDER RUNGE-KUTTA 



MATRIX INVERSION 

SOLVE LINEAR EQUATIONS - ONE R H* SIDE 
MATRIX DECOMPOSITION Ax=b 
FORWARD AMD BACK SOLUTION Ax=b 
MATRIC DECOMPOSITION BY GROUTS METHOD 
TRIANGULAR MATRIX PRODUCT LU=A 
INVERT QUASI LOWER TRIANGULAR MATRIX 
INVERT QUASI UPPER TRIANGULAR MATRIX 
TRIANGULAR MATRIX PRODUCT UL=A 
procedure ATOLIU 
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LE 11 procedure LIUSOL 

LE 12 procedure INVL1 

LE 13 pro oodu fo INVTI 

LE 14 procedure UL1TQA 

LE 15 SOLVE" TRIBIAGQNAL LINEAR EQUATIONS 

LE 16 DECOMPOSE SYMMETRIC POSITIVE DEFINITE MATRIX 

LE 17 INVERT LOWER TRIANGULAR MATRIX IN SITU 

LE 18 MATRIX MULTIPLICATION 

LE 19 DECOMPOSE BANDMATRIX INTO LOWER AND UPPER TRIANGLES 

LE 20 FORWARD SOLUTION AND BACK SUBSTITUTION FOR BAND EQUATION Ax=b 

LE 21 SOLVE BA!H) EQUATIONS Ax=b 

LE 22 REARRANGE PERSfJTED TRIANGULAR MATRIX 

LE 23 p^ocsdure I30VEL 

LE 24 DECOMPOSE SYMMETRIC POSITIVE DEFINITE MATRIX 

LE 25 FORWARD AND BACK SOLUTION Ax=b, SYMMETRIC A 

LE 26 DECOMPOSE BAND MATRIX WITH PIVOTING 

LE 27 FORWARD AND BAGS SOLUTION Ax=b , BANDMATRICES 

LE 23 DECOMPOSE SYMMETRIC POSITIVE DEFINITE MATRIX 

LE 29 FORWARD AND BACKWARD SOLUTION OF Ax=b, SYMMETRIC 

LL 01 EIGENVALUES AND EIGENVECTORS OF SYMMETRIC MATRIX 

LL 02 SOLUTION OF TEE E IGENPROBLEM (A-Ab)x=0 

LZ 01 GENERATE UNSYKMETRIC T3S MATRIX 

MC 01 LEAST SQUARES POLYNOMIAL FIT 

MD 01 LOCATE MINIMUM OF FUNCTION f (x) 

MS 01 CONVERT SEXAGESIMAL ANGLES TO RADIANS 

MG 02 CONVERT RADIAN ANGULAR MEASURE TO SEXAGESIMAL UNITS 

MG 03 CONVERT CENTESIMAL ANGLES TO RADIANS 

MG 04 CONVERT RADIAN ANGULAR MEASURE TO CENTESIMAL UNITS 

KG 05 COMPUTE DISTANCE AND GRID BEARING 

MG 06 COMPUTE POLAR CO-ORDINATES OF POINT X,Y 

MG 07 COMPUTE POINT OF INTERSECTION FROM ANGLES 

MG 08 COMPUTE POINT OF INTERSECTION FROM BEARINGS 

MG 09 THREE POINT RESECTION FY ANGLES 

MG 10 MULTIPLE POINT RESECTION FROM ANGLES 

MG 11 THREE POINT RESECTION FROM DISTANCES 

MG 12 MULTIPLE POINT RESECTION USING DISTANCES 

MG 13 THE INACCESSIBLE BASE PROBLEM 

KH 01 COKPUTS FOURIER COEFFICIENTS 

MH 02 SUM FOURIER SERIES 



MR 01 LOCATE ROOT OF f (x)=9 



KS 01 GENERAL SUM SERIES 

MS 02 ADAPTIVE SIMPSON INTEGRATION 

MS 03 EVALUATE DEFINITE INTEGRAL 



MT 01 LAGRANGE INTERPOLATION 

MT 02 AITKEN UNEQUAL INTERVAL INTERPOLATION 

MT 03 AITKEN EQUAL INTERPOLATION 

MT 04 AITKEN Nth ORDER INTERPOLATION 

MT 05 ESTIMATION OF DERIVATIVE OF GIVEN FUNCTION - UNEQUAL INTERVALS 

06 ESTIMATION OF DERIVATIVE OF GIVEN FUNCTION EQUAL INTERVALS 
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NC 01 REARRANGE ELEMENTS OF VECTOR 

NC 02 PERMUTE ROWS OR COLUMNS OF MATRIX 

NC 03 PERMUTATION OF ELEMENTS OF A VECTOR 

NC 04 PERMUTE RG17S OR COLUMNS OF MATRIX 

NC 05 PRE OR POST MULTIPLY MATRIX BY PERMUTED IDENTITY MATRIX 

NC 06 I N VERS EPERMUAT I ON OF INTEGER VECTOR 

OP 01 CRITICAL PATH SCHEDULING 

SS 01 REAL RANDOM NUMBER GENERATOR 

SS 02 REAL RANDOM NUMBER GENERATOR 

SS 03 RANDOM NUMBER GENERATOR - POISSON DISTRIBUTION 

SS 04 RANDOM NUMBER GENERATOR - STANDARD NORMAL DISTRIBUTION 

05 GENERATE LARGE INTEGER RANDOM NUMBERS 

SS 06 GENERATE LARGE INTEGER RANDOM NUMBERS 

SS 07 GENERATE SMALL INTEGER RANDOM NUMBERS 

SS 03 GENERATE SMALL INTEGER RANDOM NUMBERS 

ZM 01 OUTPUT BINARY VALUE 

ZM 02 OUTPUT OCTAL VALUE 

ZM 03 SENSE NUMBER GENERATOR KEY 

ZM 04 INPUT NEWLINE STRINGS 



ZP 02 PRINT MATRIX 

ZP 03 CONTROL FLEXOWRITER PAGE PRINTING 



TITLE 



REFERENCES 



PERIODIC PAYMENT AS FN OF UNPAID PRINC 
SORTING- MATHSORT 

PARTITION-USED BY QUICKSORT ACM64 

QUICKSORT • 

FIND 

SORTING PROCEDURES 

TREESORT- SUPERSEDED BY TREESORT 1 ACM143 
TREE SORT 1 
TREESORT 2 

SHUTTLE SORT-EXCHANGING PAIRS 

S HELL SORT- SHELLS METHOD 

STRINGSORT 

RADIX EXCHANGE ETC 

SORTING BY DISTRIBUTION COUNTING 

HEAP SORT 

TREESORT 3 

SEARCH IN A LIST 

INSERTION IN A LIST 

DELETION FROM A LIST 

SORTING WITH MIN STORAGE 



COMPUTE CODE STRING TO MOVE PLOTTING PEN 
PROCEDURES FOR GRAPH PLOTTING 



CA 

ACM45 4-61 9-63 
DH 

ACM23 11-60 5-61 
ACM63 7-61 8-62 8-63 
ACM64 7-61 8-62 8-63 
ACM65 7-61 8-62 8-63 
ACM76 1-62 6-62 
ACM113 8-62 
ACM143 12-62 
ACM144 12-62 

ACM175 6-63 10-63 12-63 5-64 

ACM201 8-63 6-64 

ACM207 10-63 10-64 

ACM 5-63 P209 

ORN13 148 

ACM 232 6-64 

ACM245 12-64 

JACM 1962(23) 

JACM 1962(23) 

JACM 1962(24) 

JACM 1962(27) 

DP 

ACM162 4-63 8-63 8-64 

NPL 0003 



ADD ITEM TO CHAIN LINKED LIST 
REMOVE ITEM FROM CHAIN-LINKED LIST 
ENLARGEMENT OF A GROUP 
LOCATE VECTOR IN LEX I COG ORDERED LIST 



DZ 

ACM100 6-62 
ACM101 6-62 
ACM136 11-62 
ACM151 2-63 



BESSEL FN 1-SERIES EXPANSION 
BESSEL FN 1 -ASYMPTOTIC EXPANSION 
BESSEL FN FOR SET OF INTEGER ORDERS 
RICCATI-BESSEL FNS OF 1ST AND 2ND KIND 
BESSEL FUNCTIONS BY RECURSION 
EVALUATE BER OR BEI FN 
Q BESSEL FUNCTIONS 

BESSEL FUNCTIONS-CHEBYSHEV-CLENSHAV7 
BESSEL FCNS 1 ST KIND 
Q BESSEL FCN 



FB 

ACM5 4-60 
ACM6 4-60 
ACM21 11-60 
ACM22 11-60 
ACM44 4-61 

ACM57 4-61 7-62 8-62 
ACM214 11-63 6-64 
AIS13 

ACM236 8-64 
ACM228 5-64 



COMPLEX EXPONENTIAL INTEGRAL 
EXPONENTIAL OF COMPLEX NUMBER 
LOG OF COMPLEX NUMBER 
NTH ROOTS OF A COMPLEX NUMBER 
COMPLEX NUMBER TO A REAL POWER 
COMPLEX DIVISION 
COMPLEX ARITHMETIC 
COMPLEX POWER 

HYPERGEOMETRIC FN-COMPLEX PARAMETERS 
CONFLUENT HYPERGEOM FN-COMPLEX PARAMS 
FREQUENCY RESPONSE 

ARSENAL OF COMPLEX ALGOL PROCEDURES 
COMPLEX SQUARE ROOT 
ZERO OF COMPLEX FUNCTION 
MULTIPLY COMPLEX VALUES 
COMPLEX PART OF POLYNOMIAL 
LOG OF COMPLEX NUMBER 



FC 

ACM« 7-60 2-61 

ACM43 4-61 6-62 

ACM48 4-61 6-62 7-62 8-64 

ACM53 4-61 7-61 

ACM106 7-62 11-62 

ECM116 8-62 

ACM186 7-63 

ACM190 7-63 

ACM191 7-63 4-64 

ACM102 7-63 4-64 

AIS16 

BIT 2 4 P237 
MCA XI AP216 
MCA XI AP217 
MCA XI AP218 
MCA X1 AP219 
ACM243 11-64 



GENERAL ORDER ARITHMETIC 
EXPONENTIATION OF SERIES -COEFFS 
EXPONENTIATION OF SERIES-COEFFS 



COMPLETE ELLIPTIC INTEGRAL 1ST KIND 

COMPLETE ELLIPTIC INTEGRAL 2ND KIND 

INCOMPLETE ELLIPTIC INTEGRALS 

COMPLETE ELLIPTIC INTEGRAL 

COMPLETE ELLIPTIC INTEGRALS 

PERIOD OF ELLIPTIC FUNCTION 

ARSENAL ELLIPTIC INTEGRALS AND FUNCTIONS 



CHEBYSHEV POLYNOMIAL BY RECURSION 
KERMITE POLYNOMIAL BY RECURSION 
LAGUERRS POLYNOMIAL BY RECURSION 
LEGENDRE POLYNOMIAL BY RECURSION 
POLYNML ORIGIN AND AXIS TRANSFORMER 
EVALUATE TABLE OF CHEBYSHEV POLYNMLS 
REDUCE DEGREE OF APPROXG POL-TELESCOPE 1 
REDUCE DEGREE OF APPROXG POL-TELESCOPE 2 
COEFFS OF QUOTIENT OF 2 POWER SERIES 
REVERSION OF SERIES 
CHEBYSHEV COEFFS 
CHEBYSHEV POLY COEFF 
POISSON-CHARLIER FOLY 



COMPLEX EXPONENTIAL INTEGRAL 
GAMMA FUNCTION 
GAMMA FUNCTION 

ASSOC LEGENDRE FNS 1ST KIND-REAL OR IMAG 

SPHERICAL NEUMANN FUNCTION 

GAMMA FN-RANGE 1 TO 2 

SET ASSOC LEGENDRE POLYNMLS 2ND KIND 

RECIPROCAL GAMMA FN-REAL ARGUMENT 

JACOBI SYMBOL BY QUADRATIC RECIPROCITY 

REAL ERROR FUNCTION 

HANKSL FUNCTIOII 1ST KIND 

PSIF-LGG DERIV OF FACTORIAL FUNCTION 

MODIFIED KANKEL FUNCTION 

INCOMPLETE BETA RATIO 

ERROR FUNCTION-LARGE X 

COMPLEMENTARY ERROR FN-LARGE X 

FIYPSRGEOMETRIC FIT-COMPLEX PARAMETERS 

CONFLUENT EYPFRGSOM FN-COMPLEX PARAMS 

CALENDAR DATE TO AND FROM JULIAN DAY 

GAUSS FUNCTION 

FRESNEL INTEGRALS 

CALCULATION OF EASTER 

COMPUTING OF ROMAN FUNCTION 

CALCN OF NORMAL DISTRIBUTION FUNCTION 

COMPLEMENTARY FRESNEL INTEGRALS 

COMPUTATION OF THE FERMI FUNCTION 

CONTD FRCTN EXPSN FOR BINML QUADTC SURDS 

EXPONENTIAL INTEGRAL BY EPSILON ALGTHM 

GRAMMA FCN 

INCOMPLETE BETA FCN 

GAMMA FCN CONTROLLED ACCY 

NORMAL DIST FCN 



FD 

ACM93 6-62 10-62 
ACM134 11-62 7-63 
ACM158 3-63 7-63 9-63 

FE 

ACM55 4-61 4-63 
ACM56 4-61 

ACM73 12-61 10-62 2-63 
ACM149 12-62 4-63 
ACM165 4-63 
NUM 4-5 P396 
NUM 5-4 P296 

FP 

ACM10 6-60 4-61 
ACM11 6-60 
ACM12 6-60 
ACM13 6-60 4-61 
ACM29 11-60 
ACM36 3-61 

AC107 3-61 8-62 8-63 
AO08 3-61 8-63 
ACM131 11-62 
ACM193 7-63 12-63 
ACM 2-62 P93 
ACM227 5-64 
ACM234 7-64 

FS 

ACM13 7-60 2-61 

ACJS1 2-61 12-62 1-63 

ACM34 2-61 7-62 

ACM 7 4-61 8-63 

AO/149 4-61 

ACM54 4-61 

ACM62 7-61 12-61 

ACM80 3-62 

ACM99 6-62 11-62 

ACM123 9-62 6-63 U-63 3-64 

ACM124 9-62 

ACM147 12-62 4-63 

ACM163 4-63 9-63 

ACM179 6-63 

ACM180 6-63 

ACM181 6-63 

ACM191 7-63 

ACM192 7-63 

ACM199 8-63 

ACM209 10-63 3-64 

ACM213 10-63 

ACM4-62 P209 11-62 P556 

ALY 1 1 P25 

ALY 1 1 F33 

BIT 2 3 PI 92 

BIT 3 2 P141 

NUM 5-2 P115 

MCA MR60 

ACM 221 3-64 10-64 
ACM 222 3-64 4-64 
ACM 225 5-64 10-64 
ACM226 5-64 
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ARCCOSSIN-ARCCOS AND ARCSIN 
TAN-LAMBERTS CONTINUED FRACTION 
ARCTAN ALD 
ARCTAN ALS 

SELECTIVE SUMMATION OF FOURIER SERIES 
ELEMENTARY FCNS CONT. FRAC, 
ARCTAN(Z) 



FT 

ACM206 9-63 
ACM 2-62 P89 
MCA XI AP1 16 
MCA XI AP1 17 
CJ 6-3 P24S 
ACM229 5-64 
ACM241 9-64 



PRODUCT OF TERMS OF A FN-CONTINUED PROD 
ALGORITHMS FOR SPECIAL FUNCTIONS 1 



FZ 

MCA XI AP202 
NUM 4-5 F409 



SECOND ORDER EQN-ADAMS METHOD 
INTEGRATION OF ORDINARY DIFFERENTL EQNS 



GA 

ACM 2-62 P88 
MCA R743 



RUNGE KUTTA 

ZERSOL-ZEROS OF SOLN OF 1ST EQUNS 
KUTTA MERSON 
FREQUENCY RESPONSE 
RUNGE KUTTA 



GL 

ACM9 5/60 
ACM194 8-63 
ACM218 12-63 10-64 
AIS16 

NUM 2 2 P131 



MULTIREGION TOO GROUP PROGRAMS 
USE O? METHOD OF KERNEL FUNCTION 
SOLVE 5 FT APPRXN DIRICHLET PROBLEM 
SOLVE 9 PT APPRXN DIRICHLET PROBLEM 
SOLVE 5 PT APPRXN DIRICHLET PROBLEM 2 
INTEGRATE HEAT EQUATION 

CANONICAL COEFFICIENTS FOR DIRICHLET PRO 



GP 

AIS14 

NUM 3 3 P219 
B0001 P49 
BOO'Jl P49 
B0001 P55 
B0001 P83 
NPL0001 



REDUCTION OF ARRAY TO JACOBI FORM 
REAL SYMM MATRIX TO TRIDIAG FORM 
GRTHONORMALISING PROCEDURE 
SYMM BANDMATRIX TO TRIDIAG FORM 
ORTHONQRMALISE COLUMNS 

REDUCTION OF A MATRIX TO CODIAGONAL FORM 
TRANSFORM MATRIX SYMMETRIC TO TRIDIAG 



LA 

ACM104 7-62 
ACM122 9-62 3-64 
ACM127 10-62 
ACM183 6-63 
ACM 2-62 P91 
CJ 4 P175 
MCA XI AP210 



ANCESTOR 
SHORTEST PATH 
PATH MATRIX 



LB 

ACM96 6-62 3-63 
ACM97 6-62 
ACM141 11-62 
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SIM EQNS-CRQUT WITH PIVOTING 
T*UDIAGOMAL SIM EQNS 
TRIDIAGQNAL LINEAR EQUATIONS 
INVERT MATRIX 

SIM EQNS-CROUT WITH PIVOTING 2 

INVERSE OF FINITE SEG OF HILBERT MATRIX 

ADJUST INVERSE OF PERTURBED MATRIX 

MATRIX INVERSION-GAUSS- JORDAN 

INVERSION-SQ RT METHOD 

SIM EQUNS AND INVERSION - BY ITERATION 

SIM EQNS-GAUSS METHOD 

MATRIX INVERSION 2-GAUSS JORDAN . 

SIM EQNS-GAUSS METHOD 

CROUT WITH EQUILIBRATION AND ITERATION 

MATRIX INVERSION 

INVERSE OF SYMMETRIC MATRIX 2 

1 ROW OF INVERSE MATRIX- MONTE CARLO 

BANDSOLVE-BAND EQUATIONS 

GAUSS SEIDEL 

INVERSION-GAUSS JORDAN 

BAND MATRICES 

GAUSSIAN SOLUTION OF LINEAR EQUATIONS 

SOLVE LINEAR EQUATIONS 

INVERT 

DETERMINANT AND INVERSE 

SOLVE USING AP220 EVAL OF DETERMINANT 

INVERT USING AP220 

INVERT USING AP220 2 

CONJUGATE GRADIENT METHOD 

SOLVE USING AP224 

INVERT USING AP224 

INVERT USING AP224 2 

SOLVE USING AP228 

GAUSSIAN ELIMINATION 

SOLVE BY QVE ^RELAXATION 

SOLVE BANDMATRIX 

SOLVE DETERMINANT INVERT SYMMETRIC 
SOLVE DETERMINANT SYMMETRIC BAND 
MATRIX PERMUTATION 
MATRIX INVERSION COMP. DIV. 
CONT. GRAD» METHOD 



FORM SET OF TEST MATRICES 
FORM SET OF TEST MATIECES 



SYMMETRIC MATRICES- MATRIX SCHEME 



MATRIX DIVISION-SQUARE ROOT METHOD 
TERM BY TERM ARITHUETRIC PROCEDURES 
SYMMETRIC MATRICES -MATRIX SCHEME 



MATRIX EIGEN VALUES AND VECTORS-JACQBI 
PREPARE FOR GIVENS METHOD AP212 
GIVENS METHOD 
COMPUTE EIGENVECTOR 
TRANSFORM EIGENVECTOR 

EIGENVALUES AND VECTORS OF SYM MATRIX 
HOUSEHOLDERS METHOD FOR SYMMETRIC MATRIX 
E VALUES OF SYM TRIDIAG MAT BY BISECTION 
EVECS OF SYM TRIDIAG MAT BY INV ITERATN 
EVALUES OF SYM BAND MATRIX BY C HOLE SKI 
EIGENVALUES AB MINUS K1 -SYMMETRIC 
EIGENVALUES A MINUS KB- SYMMETRIC 



LE 

ACM16 9-60 10-60 3-61 
ACM17 9-6C 
ACM24 11-60 

ACM42 4-61 11-61 1-63 8-© 

ACM43 4-61 8-63 

ACM50 4-61 1-62 1-63 

ACM51 4-61 7-62 

ACM58 5-61 6-62 8-62 12-62 

ACM66 7-61 1-62 6-62 

ACMS2 5-62 

ACM107 7-62 1-63 8-63 
ACM120 3-62 1-63 8-63 
ACM126 10-62 
ACM135 11-62 7-64 
ACM140 11-62 8-63 
ACM150 2-63 7-63 3-64 
ACM166 4-63 9-63 
ACM195 8-63 
ACM220 12-63 6-64 
ACM 2-62 P94 
AIS11 

BIT 2 4 P256 
MCA XI AP205 
MCA XI AP206 
MCA XI AP208 
MCA XI AP221 
MCA XI AP222 
NUM5-2 P194 
NUM MR63 AP225 
MCA MR63 AP225 
MCA MR63 AP226 
MCA MR63 AP227 
MCA MR63 AP229 
B0001 P25 
B0001 P26 
BIT 3-3 P207 
NPL 0004 
NPL 0005 
ACM230 6-64 
ACM231 6-64 
ACM238 8-64 

LF 

ACM52 4-61 8-61 11-61 8-62 
ACM52 1-63 8-63 

LG 

AIS21 
LH 

ACM197 8-63 3-64 

AIS20 

AIS21 

LL 

ACM85 4-62 8-62 8-63 

MCA XI AP211 

MCA X1 AP212 

MCA XI AP213 

MCA XI AP214 

MCA XI AP215 

NUM 4-4 P357 

NUM 4-4 P363 

NUM 4-4 P371 

NUM 5-3 P277 

NPL 0007 

NPL 0006 
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DETERMINANT BY TRIANGULARISATION 
DETERMINANT-COMBINATORIAL METHOD 
REDUCTION OF MATRIX WITH FOLYNML ELEMS 
DETERMINANT OF SYMMETRIC POSITIVE MATRIX 
FORM DETERMINANT 
DETERMINANT AND SOLUTION 
DETERMINANT AND SOLUTION 
DETERMINANT AND INVERSE 
DETERMINANT OF A BAND MATRIX 
DETERMINANT OF SYMMETRIC POSITIVE MATRIX 
DETERMINANT OF SYMMETRIC POSITIVE MATX 2 
EVALUATION OF DETERMINANT 



LQ 

ACM41 4-61 9-63 3-64 
ACM159 3-63 12-63 
ACM170 4-63 8-63 7-64 
MCA MR €3 AP224 
MCA XI AP 204 
MCA XI AP207 
MCA XI AP2Q7 
MCA XI AP20S 
MCA XI AP209 
MCA XI AP220 
MCA MR63 AP228 
ACM224 4-64 12-64 



CRAM-READ SQUARE SYMM MATRIX AS VECTOR 



LR 

ACM67 7-61 6-62 



SYMMETRIC MATRICES-MATRIX SCHEME 



LZ 

AIS21 



LEAST SQUARE FIT BY ORTHOGONAL POLYNOMIALS 
CURVE FITTING WITH CONSTRAINTS 
CHEBYSHEV CURVE FIT 

SURFACE FIT -ORTHOGL POLYNOMIAL LST SQS 

LEAST SQUARES SURFACE FIT 

LEAST SQUARES SOLUTION WITH CONSTRAINTS 

ERLANG PROBABILITY FOR CURVE FITTING 

NORMAL PROBABILITY FOR CURVE FITTING 

LEAST SQUARES SOLUTION 

SPLINE CURVE 



MC 

ACM23 11-60 12-61 
ACM74 1-62 6-63 
ACM91 5-62 4-63 5-64 
ACM164 4-63 3-63 
ACM176 6-63 
ACM177 6-63 7-63 
ACM184 7-63 
ACM185 7-63 
ACM 2-62 P92 
BIT 2 2 P81 



INTERPOLATN-DIFFERENTIATN AND INTEGRATN 
MIN OF FN OF N VARABLES- STEEPEST DESCENT 
DIRECT SEARCH-MIN OF FN OF N VARIABLES 
DIFFERENCES AND DERIVATIVES 
STEEP 1-MIN BY STEEPEST DESCENT 
STEEP 2-MIN STEEPEST DESCENT 
ATIVE-USED BY ACM 203 
MAXIMUM OF A FUNCTION 



MD 

ACM77 2-62 6-62 8-63 11-63 
ACM129 11-62 9-63 
ACM178 6-63 
ACM187 7-63 
ACM203 9-63 10-63 
ACM204 9-63 
ACM205 9-63 
MCA XI AP201 



PROCEDURES FOR RANGE ARITHMETIC 
CALCN OF NORMAL DISTRIBUTION FUNCTION 
ALGOL PROCEDURES FOR RANGE ARITHMETIC 



ME 

ACM61 7-61 
ALY 1 1 P33 
SUC 15 



POSITION OF POINT WRT POLYGON 
COMPUTE CODE STRING TO MOVE PLOTTING 
PERSPECTIVE DRAWING PROGRAM-MARK 3 



PEN 



MG 

ACM112 8-62 12-62 
ACM162 4-63 8-63 
AIS4 



SUMMATION OF FOURIER SERIES 
FOURIER SERIES APPROXIMATION 
EPSILON ALGORITHM 
EPSILON ALGORITHM! 



MH 

ACM128 10-62 7-64 
ACM157 3-63 9-63 10-63 
11-63 5-64 

NUM MATH V6 
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ROQTFINDER-FAST ITERATION METHOD 

ROOTS OF POLY-MIR5TOW-HITCHCGCK METHOD 

ROOTS BY ITERATED BISECTION 

ROOTFINDER 2- ITERATION 

iiEAL ZEROS OF AN ARBITRARY FUNCTION 

ROOTFINDER 3 

ROOTS OF POLY-NEWTON BAIRSTOW ETC 

ZEROS iiEAL POLYNML BY RESULTNT PROCEDURE 

FACTORS OF POLYNOMIALS 

ROOTS OF POLYNOMIALS- INT-COEFFS 

ZEROS OF POLYNOMIALS-NEWTON MAEHLY 

BOUNCS ON A ZERO OF A POLYNOMIAL 

.iOOTS OF ARBITRARY FN-MULLERS METHOD 

ROOTFINDERS-LINEAR AND QUADRATIC 

ZERO OF COMPLEX FUNCTION 

GRAIFEE PROCESS 

ZERO OF A FUNCTION 

SMALL INTEGERS AS PRIMITIVE ROOTS 

QUOTIENT DIFFERENCE METHOD 

ZERO OF AN ARBITRARY FUNCTION 



QUADRATURE I NT- SEVERAL FNS-SAME LIMITS 
EULER SUMMATION OF SERIES 
REAL EXPONENTIAL INTEGRAL 
MYLTIPLE INTEGRATION-GAUSS 
ROMBEilG INTEGRATION 

INTERPOLATN-DIFFERENTIATN AND INTEGRATN 
SIMPSONS INTEGRATION 

FRESNEL SIN AND COS INTEGRALS-ASYMP EXPN 

FRESNEL SIN INTEGRAL 

FRESNEL COS INTEGRAL 

DEFINITE COMPLEX LINE INTEGRALS 

SIMPSONS RULE INTEGRATOR 

DEFINITE EXPONENTIAL INTEGRALS A 

DEFINITE EXPONENTIAL INTEGRALS B 

WEIGHT COEFFS FOR GAUSS QUADRATURE 

ADAPTIVE INTN BY SIMPSONS RULE 

MULTIPLE INTEGRATION 

NONRECURSIVE ADAPTIVE INTEGRATION 

ADAPTIVE AND MULTIPLE INTEGRATION 

GAUSS FUNCTION 

FRESNEL INTEGRALS 

COMPLEIvlENTARY FRESNEL INTEGRALS 

THE LIMIT OF A CONVERGING SEQUENCE 

SIMPSON INTEGRATION VARIABLE H 

INTEGRATION BY SIMPSONS RULE 

SIMPSONS RULE IN ALGOL 58 

EULER- SUMMATION OF SERIES 

TRANSFORMATION OF SLOWLY CONVERGR SERIES 

MULTIPLE INT. SIMPS. RULE 

ROMBERG INTEGRATION 

FRESNEL INTEGRAL 



MR 

ACM2 2-60 6-60 8-60 3-61 
AOS 2-60 6-60 2-61 3-61 
ACM4 3-60 3-61 
ACM15 8-60 11-60 3-61 
ACU25 11-60 3-61 
ACM26 11-60 3-61 
AOM3C 12-60 5-61 1-62 
ACM59 5-61 

ACM75 1-62 7-62 8-62 
ACM78 2-62 3-62 8-62 
ACM105 7-62 7-63 
ACM174 6-63 
ACM196 8-63 
AIS15 

MCA XI AP217 
JCM 7 PS383 
MCA XI AP200 
NUM 4-4 E338 
MUM 5-2 P96 
BIT 3-2 P205 

MS 

ACM 2-60 

ACM8 5-6C 11-63 

ACM20 10-60 2-61 4-61 

ACB02 2-61 2-63 

ACM60 6-61 3-62 5-62 7-64 

ACM77 2-62 6-62 8-63 11-63 

ACM84 4-62 7-62 8-62 11-62 

ACM88 5-62 10-63 

ACM89 5-62 10-63 

ACM90 5-62 10-63 

ACM103 6-62 

ACM103 6-62 

ACM1G8 7-62 

ACM109 7-62 

ACM125 10-62 

ACM145 12-62 4-63 5-64 

ACM146 12-62 

ACM182 6-63 4-64 

ACM195 3-63 

ACM209 10-63 

ACM213 10-63 11-64 

BIT 2 3 PI 92 

BIT 1 1 P64 

BIT 1 4 P290 

MCA XI AP203 

NUM 1 1 P60 

NUM 2 2 P13C 

NUM 4-1 P10 

ACM233 6-64 

NUM MATH V6(15) 

ACM244 11-64 
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INTERPOLATION BY CONTINUED FRACTIONS 

INTERPOLATION BY AITKEN ITERATION 

I NTERPOLATN- DIFFERENTI ATN AND INTEGRATION 

DIFFERENCE EXPRESSION COEFFICIENTS 

CONFLUENT DIVIDED DIFFERENCES 

BACKWARD DIVIDED DIFFS 

FORWARD DIVIDED DIFFS 

DIFFERENCES AND DERIVATIVES 

LAGRANGIAN INTERPOLATION 

HERMITE INTERPOLATION 

NEVILLE INTERPOLATION 

INTERPOLATION WITH RATIONAL FUNCTIONS 



ORTHGNORMALISING PROCEDURE 
INOUT-ALGEBRA OF SETS- SPECIAL SUM 
CALENDAR DATE TO AND FROM JULIAN DAY 
SHANKS 

QRTHONQRMALISE COLUMNS 
CALCULATION OF EASTER 
SUM 

INNERPRODUCT 

FACTOR 

.REMAINDER 

SOLVE LICHTENSTEIN-GERSHOGORIN EQUATION 
CONFORMAL MAPPING OF NEARBY CIRCULAR REG 
CONVERGING FACTOR FOR CONTINUED FRACTION 
DSPLAY EPSILON ALGORITHM 
CALCULATE EPSILON ALGORITHM ARRAY 
DISPLAY SINGULAR EPSILON ALGORITHM 
CALC SINGULAR EPSILON ALGORITHM ARRAY 



BINOMIAL COEFFICIENTS BY RECURSION 
FACTORIAL 

PERMUTATIONS OF INTEGERS TO N 
COMPOSITION GENERATOR 

PERMUTE FIRST N COMPONENTS OF VECTOR X 
PERMUTATION GENERATOR-LEXICOGRAHIC ORDER 
COMBINATION GENERATOR 
PERMUTATION IN LEXICOGRAPHICAL ORDER 
PERMUTE FIRST N COMPONENTS OF VECTOR X 
PERMUTATION GENERATOR 

LOCATE VECTOR IN LEXICOG ORDERED LIST 
NEXCOM-COMBS OF ELEMNTS IN COLUMN VECTOR 
COMBINATION IN LEXICOGRAPHICAL ORDER 
COMBINATION IN ANY OuDER 
CALC COMBS OF M THINGS N AT A TIME 
COMBS FROM 1 TO N AT A TIME OF M THINGS 
PERMUTATIONS IN LEXICOGRAPHICAL ORDER 
RANDOM PERMUTATION 
PERMS WITH REPETITION 

SIEVE-PRIME NUMBERS UP TO N 
PRIME TWINS 

GREATEST COMMON DIVISOR 



GENERATE PARTITIONS IN PART COUNT FORM 
GENERATE PARTITIONS WITH CONSTRAINTS 



MT 

ACM 18 9-60 3-62 

ACM70 11-61 7-62 

ACM77 2-62 6-62 8-63 11-63 

ACM79 2-62 3-63 

ACM167 4-63 9-63 

ACM168 4-63 9-63 

ACM169 4-63 9-63 

ACM187 7-63 

ACM21G 10-63 10-63 P619 

ACM211 10-63 

AIS17 

NUM 3-3 F302 

MZ 

ACM127 10-62 
ACM156 3-63 8-63 
ACM199 8-63 
ACM215 11-63 
ACM 2-62 P91 

ACM 4-62 P209 11-62 P556 

MCA XI AP119 

MCA XI API 20 

MCA XI AP125 

MCA XI API 26 

B0001 P242 

B0001 P256 

NUM 5-4 P341 

BIT 3-3 PI 83 

BIT 3-3 PI 84 

BIT 3-3 PI 91 

BIT 3-3 PI 92 

NC 

ACM19 10-60 6-62 8-62 
ACHS3 2-61 

ACM71 11-61 4-62 8-62 
ACM72 11-61 3-62 
ACM86 4-62 8-62 
ACM87 4-62 8-62 10-62 
ACM94 6-62 11-62 12-62 
ACM102 6-62 

ACM115 3-62 10-62 12-62 
ACM130 11-62 
ACM151 2-63 
ACM152 2-63 7-63 
ACM154 3-63 8-63 
ACM155 3-63 8-63 
ACM160 4-63 8-63 10-63 
ACM161 4-63 3-63 10-63 
ACM202 9-63 
ACM235 7-64 
ACM242 10-64 
NP 

ACM35 3-61 4-62 8-62 
ACM223 4-64 
ACM237 3-64 12-64 

NR 

ACM95 6-62 
ACM114 8-62 
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EUCLIDEAN ALGORITHM 
AUGMENTATION OF X BY Y 
CHAIN TRACING 

JACOBI SYMBOL BY QUADRATIC RECIPROCITY 
MAGIC SQUARE-EVEN ORDER 
MAGIC SQUARE-ODD ORDER 
DIOPHANTINE EQUATION 
TERM OF MAGIC SQUARE 

CALENDAR DATE TO AND FROM JULIAN DAY 
CALCULATION OF EASTER 



ASSIGNMENT 

ECONOMISING A SEQUENCE 1 
ECONOMISING A SEQUENCE 2 
OPTIMAL CLASSIFICATION OF OBJECTS 
SHORTEST PATH 

INTEGER SOLNS OF LINEAR PROG PROBLEMS 
DISCRETE CONVOLUTION 



CRITICAL PATH SCHEDULING 
EVALUATION OF A PERT NETWORK 
MINIMUM EXCESS COST CURVE 
TOPOLOGICAL ORDERING FOR PERT NETWORKS 
A PERT PROGRAM IN ALGOL 60 



CARBON DIOXIDE AND AIR FUNCTIONS 
THERMODYNAMIC PROPS OF STEAM AND WATER 



MULTIREGION TOO GROUP PROGRAMS 



QUANTUM MECH INTS- SLATER TYPE ORBITALS 

MOLECULAR ORBITAL CALCULATION 

QUANTUM MECH INTS- SLATER TYPE INTEGRALS 



CHI SQUARED 

STATISTICAL ANALYSIS-CARD INPUT 



SMOOTHING BY GRAMS FORMULAS- 1 
SMOOTHING BY GRAMS FORMULAS- 2 
SMOOTH-FOURTH ORDER SMOOTHING 



CORRELATION COEFFS T7ITH MATRIX MULT 



FREQUENCY DISTRIBUTION 
MEAN STANDARD DEVIATION AND VARIANCE 
MEAN STANDARD DEV VAR AND DEG OF FREEDOM 
CALCN OF NORMAL DISTRIBUTION FUNCTION 



TRIANGULAR REGRESSION 
MULTIPLE REGRESSION- MARK 1 



NZ 

ACM7 4-60 
ACM68 8-61 11-61 
ACM69 9-61 
ACM99 6-62 11-62 
ACM117 8-62 1-63 3-63 
ACM118 8-62 12-62 1-63 
ACM139 11-62 
ACM148 12-62 4-63 
ACM199 8-63 

ACM 4-62 P209 11-62 P556 
OL 

ACM27 11-60 10-63 12-63 

ACM81 3-62 

ACM 82 3-62 

ACM83 3-62 

ACM97 6-62 

ACM153 2-63 8-63 

ACM208 10-63 

OP 

ACM40 3-61 9-61 10-62 6-64 
ACM119 8-62 
ACM217 12-63 
ACM219 21-63 
MCA MR56 

PH 

AIS12 
AIS19 

PN 

AIS14 
PQ 

ACM110 7-62 

ACM111 7-62 

ACM132 11-62 

SA 

AIS8 
AIS22 

SB 

ACM188 7-63 

ACM189 7-63 

ACM216 11-63 

SC 

ACN39 3-61 
SM 

ACM212 10-63 

AIS6 

AIS7 

ALY 1 1 933 

SR 

ACM142 12-62 
AIS18 
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POS NORMAL DEVIATE FROM RANDOM NUMBER 

GENERATE RANDOM NUMBER 

ERLANG PROBABILITY FOR CURVE FITTING 

NORMAL PROBABILITY FOR CURVE FITTING 

NORMAL RANDOM 

DISCRETE CONVOLUTION 

QUASI-RANDOM POINT SEQUENCE 



SS 

ACM121 
ACM133 

ACM184 7-63 

ACM185 7-63 

ACM200 3-63 

ACM208 10-63 

ACM247 12-64 



9-62 

11-62 12-62 3-63 



VARIANCE HOMOGENEITY TEST-BARTLETT 
MEAN STANDARD DEVIATION AND VARIANCE 
MEAN STANDARD DEV VAR AND DEG OF FREEDOM 
TWO WAY ANALYSIS OF VARIANCE 
2 TO THE K FACTORIAL ANALYSIS- YATES 



SV 

AIS5 
AIS6 
AIS7 
AIS9 
AIS10 



NESTING OF FOR STATEMENT 1 
NESTING OF FOR STATEMENT 2 
DIAGRAM 
FUSBUDGET 

COLLAPSE OWN STORAGE 
BOOLEAN OPERATOR NOT FLAG 
STORAGE ALLOCATION PROCEDURES 
REDUNDANCY CHECK 
ALGOL 60 SYNTAX CHECKER FOR THE IBM 7G900RNIJ3399 



ZA 

ACM137 11-62 
ACM138 11-62 
ACM 1-61 P55 
ACM 1-61 P64 
ACM 1-61 P64 
ACM 1-61 P71 
ACM 10-61 P444 
ACM 6-62 P340 



ASSIGN 

CALENDAR DATE TO AND FROM JULIAN DAY 
CALCULATION OF EASTER 

ALGOL INFORMATION SHEETS-THEIR PURPOSES 

ALGOL ON PEGASUS 

ALGOL ON DEUCE 

FREE FIELD READ 

COORDS FO ELLIPSOID 

GRAYCODE 



ZZ 

ACM173 6-63 10-63 

ACM199 8-63 11-64 

ACM 4-62 P209 11-62 P556 

AIS1 

AIS2 

AIS3 

ACM239 8-64 
ACM240 9-64 
ACM246 12-64 
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LOCATE ELEMENT IN LIST 

procedure search (p) lowerbound : (£) upperbound : (u) relation : (b) ; 
integer p,£,u; boolean b; 

comment HUCC LIBRARY PROCEDURE DH02: 
AUTHOR: J. BOOTHROYD 

A Jensen procedure which searches an ordered array a 

between the inclusive limits a[l] and a[u] for 

an element of given value x, say. For arrays 

sorted in ascending order the relations 

b = (x 4 a[p]) and b = (x < a[p]) yield exit 

values of I and u satisfying a[£] < x 4 a[u] 

and a[&] 4 x < a[u] respectively. 

To locate an element of value y in the third column 
of a matrix A [ 1 : 10 , 1 : 5 ] for example, use the call 

JU-1; u:«10; 

s e a r ch ( p , il , u , y <A [ p , 3 ] ) 
and to find an element of value y in the second 
row of the same matrix perform 

£:=1; u:=5 : 

search(p,£,u,y 4 A[2,p]) 
Note from these examples that £ and u specify 
the extent of the search in the call as well as serving 
as result variables on exit from the procedure. 
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SORT REAL NUMBERS INTO ASCENDING ORDER 

procedure shellsort (a, n) ; value n; real array a; integer n; 
comment HUCC LIBRARY PROCEDURE DHOl: 
AUTHOR J, BOOTHROYD : 

Elements a[l] through a[n] of real array a[l:n] are 
rearranged in ascending order. The method is that of 
D.A. SHELL (A High Speed Sorting Procedure, COMM. A. CM. 2 
(1959) 30-32) with subsequences chosen such that m lf the first 
value of the partition parameter, is given by 

k . „ rt k k+l 
mi = 2-1 for 2 4 n < 2 

To implement Shell's original choice of m x » [n/2] change the 
first statement of the procedure to m: = n. Note that shellsort 
specifies the array as type real . It may thus not be used to 
reorder elements of integer arrays unless the specification of 
a is changed to integer array a and in this case the local 
working variable w should preferably be declared as integer w 



SORT REAL NUMBERS IN ASCENDING ORDER 



procedure exsort(a,n) ; value n; integer n; array a; 
comment HUCC LIBRARY PROCEDURE DH03: 
AUTHOR J. BOOTHROYD s 

A procedure which sorts the elements a[l] through a[n] 
of real array a [1 :n] into ascending order. The method 
used is that of exchanging element pairs. After the 
firs t pass , which ranges over all element positions 
a[ 1] to a[n] , each pass is restricted to the range 
bounded by the first and last exchanges which occurred 
on the previous pass . This significantly improves the 
efficiency of the method. 

The procedure is simple and may be easily modified 
to sort numbers into decending order by changing the 
Boolean expression a [ i+1 ] < a[i] to a[i+l] ^ a[i] . 
A version of exsort (under the identifier jensort) 
suitable for sorting rows or columns of an m x n matrix 
is included in the library as DH04 ; 



DH04 



SORT ROWS OR COLUMNS OF MATRIX 

procedure jensort (ai ,aiplusl,i,n) ; value n; real ai ,aiplusl; 
integer i ,n; 

comment HUCC LIBRARY PROCEDURE DH04: 
AUTHOR J. BOOTHROYD : 

A Jensen modification of exsort (DH03) which permits the 
sorting into ascending order of rows or columns of arrays 
of any number of dimensions. To sort the 3rd row of an 
array declared as A[l:n,l:m] use the call 
jensort(A[3,i] ,A[3,i+l] ,i,m) . To sort each column of 
B [1:10,1: 20] use the construction 

for p: = 1 step 1 until 20 do jensort (B [i ,p] ,B[i+l,p] ,i,10) . 
These examples , and a study of the differences between 
DH03 and DH04 should make clear the significance of the 
parameters ai and aiplusl; 



procedure shell8ort(a,n) ; cessment ACM 201; value n; real array a; integer n; 
begin integer i t j,k t m; real w; switch s:=endj; 
for i:=1 step i until n do m:=2*i-1; 
for m:=m div 2 while m £ do 
begin k:=n-m; 

for j:=1 step 1 until k do 

begin for i:=j step -m until 1 do 

begin if a[i+m] > a[i} then goto endj; 

w:=a[i]; a[i]:=a[i-Hn]; a[i+m]:=w 

end i; 

endj: end j 

end m 

end shellsort; 
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procedure search(p) lowerbound: (l)upper bound: <u)relation:(b) ; 
integer p 9 l 9 u; boolean b; 

comment a Jensen procedure which searches an ordered array a between the inclusive limits a[ll and a[u] for 
an element of given value, say x. For arrays sorted into ascending order the relations 
b=(xCa[p)) and b=(x<aCp3 produce exit values of 1 and u satisfying atl] < x < a[u] and 
a[l]~< x < a[u] respectively; 
begin for pT=(u+l) div 2 while u > 1 do if b then u:=p-l else l:=p+1 ; 

" u:=u-H ; 1:=1-1 
end search; 
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procedure exsort(a,n) ; value n; integer ft; array a; 
begin integer i,k,r; switch s:=L1,L2; real temp; 
Li: k:=2; 

L2: r:=n-l ; n:=0; 

for i:=k-l step 1 until r do 
if a[i+l]<a[i] then 

be S in temp:=aCil; ati]:=a[i+1]; a[i+1 ]:=temp; 
if n=0 then k:=i; n;=i 

end ; 

if n#0 then goto if k±1 then L2 else LI 
end exsort; 



procedure jensort(ai, aiplusi t i,n) ; value n; real ai f aiplust; integer 
begin integer k,r; switch s:=&1,L2; real temp; 
Li: k:=2; J 
L2: r:=n-i; n:=0; 

for ij=k-l step 1 until r do 
if aiplusi < ai then 

***giP temp:=ai; ai:=aiplusl; aiplusi t=temp; 
if n=0 then k:=i; n:=d 

end ; 

if n * then goto if kil then L2 else LI 

end jensort ; 



i 
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OUTPUT GRAPH OF VECTOR ELEMENTS 

procedure Graphic(a,m,n) ; 
value ra,n; array a; integer m,n; 
comment HUCC LIBRARY PROCEDURE DPOl: 
AUTHOR: W. G. WARNE: 

Plots the values of the elements of the nth order 
vector a on an (nrfl)Kn grid. Each element is plotted 
on a new line as a + character. Scale is chosen so that 
the maximum of aj occurs m spaces from the left, and the 
minimum zero spaces from the left. Three new line 
characters are output before the first line marked by a 

of ro+1 dots and two newlines are output after the last 
line marked by a row of m+1 dots. A column of vertical 
bar characters marks the baseline (left margin). 

Use is made of tabs to ensure optimum print speed. 
Tabs must be set at 12 spaces. 



procedure Graphic ( a , n , m) • 

value m,n; array a; j£££E£E> m > n > 

begin ijateg er j, k, h, r ; real max, min, x, y ; 

print ££13??; 

for j := steg 1 until m do print £.?; 
k := 0; 

for j := 1 s^teg 1 unti l n do 

k := i£ a£k] > at j] thon k else j; 
max := a[k] ; 
k := 0; 

for j := 1 step 1 until, n do 

k := if a[k] < a[ j] thgn k olee j; 
min := a[k]; 
x := m /(max- min) ; 
for h := step 1 until n do 
begin pri^t ££1??; 

k := x*(a[h]-min> ; 

r := k div 12; 

for j := 1 steg 1 until r do print ££t?? ; 
r s= k - r*12; ^ 
f or j := 1 stej> 1 until r do print ££c??; 
print £+?; 

end; 

print ££1??; 

for j := ste£ 1 until m do print £.?; 
p rint ££12?? 

end Graphic; 
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OUTPUT PLOT OF ELEMENTS OF TWO VECTORS 

procedure Bigraphic(a,b,n,m) ; 
value n,m; array a,b; integer m,n; 
comment HUCC LIBRARY PROCEDURE DP02 
AUTHORS W. G. WARNE: 

Plots values of the elements of two nth order vectors 
a, and b on an (m+1) x n grid. 

Each pair of elements is plotted on a new line, 
aj being represented by _ (underline character) and bj 
by a dot. Scale is chosen so that the maximum of aj or bj 
occurs m spaces from the left and the minimum of aj , bj 
at zero spaces, and both elements are plotted with the same 
scale and origin. 

Otherwise format is as for Graphic (qv) . 



[ — — , 
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pr oc ed u re Bigraphic(a f b f n,m) ; 
value n, m; arra y a,b; integer m,n; 

begin integer j^kajkb^hySpjl; 
real min, max,x; 
procedure position; 

begin integer i ; 

if l+sp<1 2 then 

begin for i:=l step 1 until sp do print ££s??; 
l:=l+sp; 
sp:=o 

end 

else if 1O0 then 

begin for i:=(l+sp)div 12 steg -1 until 1 do print ££t??; 
l:=l+sp-(l+sp) div 12 * 12; sp:=0; 
for i;=1 step 1 until 1 do print ££s?? 

<~ mju-> csbmqCs awnracs mama m 

end else 

begin for i:=12-l step -1 until 1 do print ££s??; 
° sp := sp-12+1; 1:=0; 
position 

end 

end position; 
jprint ££12??; 

for j:=0 step 1 until m do print £ # ?; 

— if aC03>b[0] then a[0] else b[0]; 
for j:=1 step 1 u ntil n do 

max:= jlf a[ j]>max or b[j]>max then (if a[j]>b[ j] then 

atj] else b[j]) else max; 
min:= i^f a[0]<b[o] then a[0] else b[0]; 
for j:=l step 1 until n do 

min: = a[j]<min or b[ jT<min then (if a[ jl<b[j] then 
f lse b C j]) else min; ~ a "~~ 
x:=m/TSax~min) ; 
for h:=0 step 1 until n do 
begin print ££!??: — — 

ka:=x*(a[h] min) ; 
kb:=x*(b[h]-min> ; 

sp:=0; . 
!:=0j I f 
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end; 



end Bigraphic; 



for j:=0 step 1 until m do 
begin if ka= j then 

begin position* 
print £? 

end; 
i^f kb = j then 
begin position; 

print £.7 

end 

else sp:=sp+1 ; 

end; 



arint ££1??; 

for j:=0 step 1 u ntil m do pri nt £,?; 
print ££!£?? 



DP02 



cont. 



) 



DZ01 

CHARACTER PACK AND UNPACK 

procedure pack(x) in position: (n)error exit: (label) ; 
comment global integer MODE equal "to the number of bits in x 

must be declared and assigned before procedure is called; 
value x,n; 
integer x,n; 
label label; 

HUCC LIBRARY PROCEDURE DZOi: 
AUTHOR: P. W, FORD: 

y 
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procedure unpack (x) from position: (n) ; 
comment as for pack procedure; 
value n; 
integer x,n; 

HUCC LIBRARY PROCEDURE DZ02: 
AUTHOR: P. W. FORD: 

This pair of procedures packs and unpacks sets of words containing 
MODE binary digits into the available ALGOL free store. In any one 
program, MODE must be constant and less than 39. The packing is dense, 
that is, each location in store is completely filled. If more than 
MODE binary digits are presented to be packed, the superfluous bits are 
ignored, n must not be less than 1. Storing starts at location 316 
and works up the store until the program's data space is reached. If 
current data space is about to be overwritten, the procedure pack outputs 
on the typewriter "STORE FULL" followed by the machine location about to 
be overwritten and the number (n) of the word which would be placed in 
the location. The procedure then exits through the label without over- 
writing the store. Procedure unpack does not check the upper bound 
since, although rubbish may be extracted, the program itself cannot be 
destroyed by this procedure. 

The machine is not aware of any packed information, so that- 
subsequent data space allocation by program could overwrite some of the 
packed store. Thus the rule is , do not declare any variables , arrays, etc. 
after procedure pack is called, at least until the packed data has been processec 
To pack and unpack one word takes of the order of two milliseconds . 

ALGOL Tape 1 is destroyed by these procedures , but the own code system 
is unaffected. 
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procedure pack(x)in position: (n) error exit: (label) ; 

comment global MODE equal to the number of bits in x must be 
declared and assigned before procedure is called; 
value x,n; 
integer x,n; 
label label; 
begin integer a,b,c,divcount, shift, count ; 
arra y A[ 1 : 1 ] ; 

switch s : terror, end, unset; 

a:=address(A)-1 ; 

b:=3l€; 

c:=39; 

elliott(%2,0,0,2,7,n) ; 

elliott (3, 0„ MODE, 0,5,2, n) ; 

elliott(5,7,o,0,0,0,0); 

elliott ( 1, 6, count, 0,5,6, c) ; 

elliott (0, 4, b,0 ? 2,0,divcount) ; 

elliott (0, 7, a, 0,4,1 : error) ; 

elliott(3,0,divcount,0,0,5,b) ; 

elliott(5,2,c,0,5, 7,0) ; 

elliott (?, 7, count, 0,0, 7, c) ; 

elliott(2,0,shift,0^6,7,divcount) ; 

elliott (3, 0,0," ,6,7, shift) ; 

elliott(5, 0,0,0, 1,%x); 

elliott(6, 7, MODE, 0,5,0,0) ; 

elliott(3,0,x,0 f 6,7,shift); 

olliott(5,4,0,0,6,7,divcount); 

elliott(2,0,0,0,0,2, shift) ; 

elliott(0, 7, MODE, 0,4,1, unset) ; 

elliott (5, 4, 39,0, 2, 0,x) j 

elliott (0, 2. divcount 5 0,0,0,0); 

elliott(l,0,x, 1,2,0,0) ; 
unset: elliott(4,3 f end,0,4,0,end) ; 

error: print punch ( 3), £ 

STORE FULL? , s ame 1 i ne , di vc ount , n ; 

goto label; 
©nd: end pack; 



procedure unpack(x) from position: (n) ; 
value n; 
integer n,x; 

be ^ in integ er a, b-c,divcount, shift, count, diff, temp; 
s ^itch ss:=rt,pu, collate, next; 
b:=316; 
c:=39; 

elliott (« ,2, 0,0, 2,1, temp); 

elliott(3,0,c,0,0,5,MODE) ; 

elliott ( 1 ,0, temp, 1,5,1,0) ; 

elliott (2,0, a, 0,3,0, n) ; 

elliott(5, 2, MODE, 0,5,7,0) ; 

elliott(l,6,count,0,0,0,0) ; 

elliott (5,6, c, 0,0,4, b) ; 

elliott(2,0,divcount,0,0 ? 5,b) ; 

elliott (5, 2,c,n,5,7,0) ; 

elliott(0,7,count,0,4,2,pu) ; 

elliott(2, 0, diff ,0,^,7, c) ; 

elliott (2 , 0, shift ,0,6, 7, divcount) ; 

elliott( 3, 0,0, 0,2,0, temp) ; 

elliott(0,2,diff ,0,0, 7, MODE) ; 

elliott(4, l.rt, 0,3,0, temp); 

elliott(5,0, 39,0,6, 7, divcount) ; 

elliott (2, 7, 81 91, 0,6, 7, diff); 

elliott(5, 4,?, 0,4,0, collate) ; 
rt: elliott (3,0, temp, 0,6, 7, shift); 

elliott(5, 1,0,0,4,0, collate) ; 
P u: elliott(0,0, divcount, 1,2,7.8191); 

collate: elliott(2,3,a,0,4,3,next); 
next: x:=a 
end unpack; 
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Double Length Floating Point Arithmetic Package . 
Author : W.G. Warne 

Procedures DZ03 to DZll provide facilities for performing real arithmetic and 
yield results which have twice the accuracy of the standard Algol and arithmetic 
operations. 

Within this scheme a double length floating point number is represented by a 
number pair each of whose members is of type integer . For the number pair 
(al, a2) , for example, the floating point number a*2 is distributed within 
al, a2 as follows 

a is held as bits 39 to 1 of al and bits 38 to 19 of a2. 

b is held as bits 10 to 1 of a2. 

bits 39 and 18 to 11 of a2 are unused. 
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With this representation the double length floating number range is 



-2 511 < x < - 2" 512 ,0,2 



-512 .511 fc4 . . 

< x < 2 with a mantissa precision 



of 57 binary digits, roughly 17 decimal digits. 



No provision is made for printing the decimal equivalent of a double length 
floating point number. The procedures are intended for special programming 
occasions when it is essential to minimise round-off errors. These procedures 
will be available only to those who convice either the of f icer-in-charge or 
the second of f icer-in-charge that their particular problem requires such precision. 
This restriction is imposed for the following reasons :- 

(a) TIME An average time for all procedures is 1 to 1,5 milliseconds. 

(b) STORAGE Total storage for the package is 659 locations. 



ERROR INDICATIONS 

Where the result of an operation exceeds the stated range, and also when the 

single length result of a package operation exceeds normal single length range, 

an error indication is given and the program is stopped. A.n error indication 

—256 

is also given if single length underflow is encountered (0<|x| <2 ) • 
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procedure expand (s) single length real 


to double 


i (dl,d2); 


value s ; real s ; integer dl,d2; 
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procedure intexp (i) becomes : (dl t d2) ; 


value i ; integer i, dl 
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procedure contract (dl,d2) double to single : (s) ; 
value dl,d2 ; integer dl,d2; real s; 
comment error is "contract error" ; 



DZ06 

procedure mpyssd (si) and :(s2) multiplied to give double :(dl,d2); 
value sl,s2; integer dl,d2; real sl,s2; 

DZ07 

procedure divdss (dl,d2) divided by: (s) gives: (r) ; 
valu e dl,d2,s ; reatl s ,r; integer dl,d2; 
comment error is "divdss error"; 

DZ08 

procedur e negd (dl,d2) has complement :(nl,n2) ; 
value dl,d2; integer dl,d2,nl,n2; 

DZ09 

procedure addddd (al,a2) plus: (bl,b2) gives: (rl,r2); 
value al,a2,bl,b2; integer al,a2,bl,b2,rl,r2; 

DZ10 | 

procedure mpyddd (al,a2) multiplied by :(bl,b2) gives :(rl,r2) ; 

valu e al,a2,bl,b2; integer al,a2 ,bl ,b2 ,rl,r2 ; 

comment error is "dfptoobig" ; 

DZ11 " "1 

procedure divddd(al,a2) divided by: (bl,b2) gives :(tl,t2); 
value al,a2,bl,b2; integ er al,a2,bl,b2,tl,t2; 

comment in the result (tl,t2) of this procedure the mantissa accuracy is 38 
binary digits, (about 12 decimal digits). 
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BESSEL FUNCTION FOR SET OF INTEGER ORDERS 

procedure BESSEL(x,n,eps ,J) ; value x ,n,eps ; real x,eps ; 
Integer n ; real array J ; 

comment HUCC LIBRARY PROCEDUPE FB01: 
SOURCE : ACM 21: 

Computes the values of the Bessel functions Jp(x) 
for real x for the set of all integer orders 0<p<n. 
These values are stored as J[p] in array J[0:n] . S 
also Comm ACM 3 (1960) p 600 and 8 (1964) p 219; 



rec; 



norm: 



exit: 



procedure BESSEL(x,n,eps,J) ; value x,n,eps; real x,eps; integer n; real array J; 

comment Based on Algorithm 21, Communications of the A.C.U. , 3, (I960) , page 600, as 
corrected and certified Vol. 8(1964) , page 219; 

begin real dist,rec0, reel ,rec2, sum, max, MAX, abx, err, z; 

integer m,k,p,q; boolean s; real array Jbar[C:n]; 

switch SS:=exit, rec, norm; 

if x=0.0 then begin J[0]:=1.G; 

for p:=1 step 1 until n do J[p]:=0.0; goto exit 

end ; 

abx:=abs(x); 

dist:= if abx > 8,0 then 5.0 *abxt(1. 0/3.0) 

else 10. G; 

k:=entier(( if abx > n then abx else n)+dist)+1 ; 
•s:= false ; MAX : =0 •*25 10 +74 ; 
recO :=rsum: =0,0; rec 1 : =1 . ; 

z:=MAX*abs(x/k) ; 

for p:=k step -1 until 1 do 

begin if p>n+1 then q:=n else q:=p-1; 
J[q3 : =rec2 : =2 . 0*p/x*ree1 -recO ; 
if abs(rec2)>z then 

begin reel : =r eel /z; rec2:=rec2/z; sum:=sum/z; 

for m:=n step -1 until p-1 do J[m]:=J[m]/z 

end; 

if p=1 then sum:=sum+rec2 

else if (p div 2*2)=fcp then sum:=sum+2.0*rec2; 
recO : =rec1 ; rec 1 : =rec2 ; 

end recursion; 
for p:=0 step 1 until n do J[p]:=J[p]/sum; 
if s then begin max:=0; 

for p:=0 step 1 until n do 

begin err:-abs(J[p]-Jbar[p]) ; 

if err>max then max:=err 
end maximum error; 
if max < eps then goto exit; 

end 

else s:= true ; 

for p:=0 step 1 until n do Jbar[p]:=J[p]; 
k:=entier(k+dist); goto rec; 
end BESSEL; 
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P, WYNN'S ARSENAL OF COMPLEX ARITHMETIC PROCEDURES 



Procedures FCOl to FC26 are taken from BIT Vol 2, No. 4, p. 237. 
Within these procedures any complex number written z in ordinary 
algebraic notation is represented by an ALGOL array of two elements, 
its real and imaginary parts. 

It should be declared thus array a[0:l] ; 
The elements z[0] and z[l] are respectively the real and imaginary 
parts. Outs ide the Wynn procedures the complex number is referred 
to as z[i] . The subscript i must be declared as an integer and 
the declaration must be valid in all blocks in which statements of 
procedures FCOl to FC26 occur. Inside the procedures i is assigned 
the value before operations on real parts, 1 before operations on 
imaginary parts. Outside the procedures i should not be used. 

Vectors , matrices and arrays with more dimensions having complex 
elements may be accommodated wi thin the general scheme. For example 
an array wi th typical element Ap,q,r might be declared as 

array A[U:ul,Jl2:u2,Jl3:u3,0:l] and 
referred to by A[p,q,r,i] 

or alternatively array A[0:l,£l:ul,£2:u2,U:u3] 
appropriately referenced by A[i,p,q,r] 

Call by name is used extensively in these procedures permitting 
the use of linear expressions of complex variables as actual parameters , 
e.g. aizi + a£Z2 + ^3^3 + ...... where a\ ,a 2 ,a 3 are real numbers and 

z l t z 2i z 3 complex scalars or complex subscripted variables. 

FCOl 

COMPLEX ASSIGNMENT 

procedure eq(one, other) ; real one, other; 

Examples: eq(z2[i] ,zl [i]) effect is "z 2 : = zi" 

eq(z3[I],zl[i]+z2[i]) "z 3 := *l+ z l 

eq(z3[i],zl[i]-z2[i]) "z 3 := zi-Z2 n 



-2- FC02 | 

MULTIPLE COMPLEX ASSIGNMENT 

procedure seqeq (third, second, f irs t) ; real third, second, first; 
Example: seqeq(z4[i] ,z3[i] ,z2 [i]+zl [i]) "z^:- z 3 := z 2 +Z!" 
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COMPLEX MULTIPLICATION 

procedure cm( res, one, other) ; real res , one, other; 
Example: cm(z3[i] ,z2 [i] ,zl[i]) "z 3 :« z 2 *z 1 " 
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COMPLEX DIVISION 

procedure cd (res, one, other) ; real res , one, other; 
Example: cd(z3[i] ,z2[i] ,zl[i]) "z 3 := z 2 /zi" 
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ASSIMILATE REAL IN COMPLEX OPERATION 

real procedure real (variable) ; real variable; 

Example: eq(z2 [i] ,zl [i]+real(a) ) "z 2 : = z^a" 



FC06 



ASSIMILATE IMAGINARY IN COMPLEX OPERATION 

real procedure imaginary (variable) ; real variable; 

Example: eq(z2[i] ,zl[i]+imaginary(a)) "z 2 := zx+j*a" 

(where, he re, j 2 = -1) 



-3- 
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COMPLEX CONJUGATE 

real procedure cxconj (it) ; real it; 
Example: eq(z2 [i] , cxconj (zl[i])) 

if zl = a + j * b then z 2 = a - j * b (j 2 - -1) 

FC08 



MODULUS OF COMPLEX NUMBER 

real procedure mod(it) ; real it; 

Example: r := mod(zl[i]) 

1 

evaluates r = | z| where |z| = (a 2 +b 2 ) 2 for z = a+j *b 

FC09 



ARGUMENT OF COMPLEX NUMBER 

real p rocedure arg(it) ; real it; 

Example: theta:= arg(zl[i]) assigns to theta the argument of z 
for z = re J the argument lies in the range -7r<0<7r 

r 
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POLAR FORM 

procedure polar forra(res ,r, theta) ; 

Example: polar form(zl [i] ,R,phi) assigns to the elements of zl the 
values rcos<J>,rsin<|>. 

FC11 



MULTIPLICATION BY IMAGINARY OPERATOR 



zi:= j*z 2 " 



Procedures FC12 to FC23 (except FC17) are all of the form 
procedure function(res ,i t) and effect the operation 
res := function(it) • 



FC12 

COMPLEX SQUARE 

procedure compsq(res ,it) ; real res, it; 

FC13 



COMPLEX RECIPROCAL 

procedure comprecip(res ,it) ; t6al res, it; 

FC14 



COMPLEX ROOT 

procedure cxsqrt (res ,it) ; real res ,it ; 
N.B. USES FC08, FC09 and FC1#. 

FC15 



COMPLEX LOGARITHM 

procedure compln(res ,it) ; real res, it; 




COMPLEX EXPONENTIAL 

procedure compexp(res ,it) ; real res ,i t ; 



FC17 



HYPERBOLIC FUNCTIONS 

procedure hyp (s inh , cosh ,y) ; value y ; real sinh,cosh,y ; 

evaluates coshy=(e^+e ^) /2 and sinhy-(e^-e ^)/2 
Special care with precision is taken in evaluating sinhy for small values 
of y. 
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COMPLEX SINE 

procedure compsin(res,it) ; real res, it; 
N.B. USES FC17. 

FC19 



COMPLEX COSINE 

procedure compcos(res,it) ; real res.it; 
N.B. USES FC17. 

FC20 



COMPLEX TANGENT 

procedure comptan(res ,it) ; real res, it; 
N.B. USES FC04 and FC17. 

FC21 



COMPLEX INVERSE SINE 

procedure cxarcsin(res ,it) ; real res ,it; 

FC22 



COMPLEX INVERSE COSINE 

procedure cxarccos(res,it) ; real res, it; 
N.B. USES FCOl, FC05, FC21. 

FC23 



COMPLEX INVERSE TANGENT 

procedure cxarctan(res ,it) ; real res.it; 

N.B. FCOl, FC04, FC05, FC12 and FC14. 

FC24 

COMPLEX POWER OF COMPLEX VARIABLE 

procedure onehochother(res, one, other) ; real res , one, other; 
Example: onehochother(z3[i] ,z2 [i] ,zl[i]) has the effect of 
z 3* = 
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TEST FOR EVEN INTEGER 

boolean -procedure even(integer) ; integer integer; 

This procedure takes the value true if 2 is a factor of integer. 
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LOGARITHM OF COMPLEX NUMBER 

procedure LOGC(a,b 5 c s d,FAIL) ; value a,b; real a,b,c,d; label FAIL 
comment HUCC LIBRARY PROCEDURE FC26: 
AUTHOR ACM 243 : 

The procedure computes the number c + di which is equal 
to the principal value of the natural logarithm of 
a + bi a i.e. such that -it < d < + ir. A label parameter 
FAIL permits exit from the pro cedure in the event that 
the real part of the result approaches minus infinity. 
Where required in the procedure the numerical values of 
12 and ln(/8) are provided: 



procedure eq(one, other) ; real one, other ; 
for i:=0,1 do one: mother; 
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procedure seqeq( third, second, first) ; real third, second, first ; 
for i :=Q, 1 do third :=second:=f irst ; 
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procedure cm( res , one , other) ; real res , one , other ; 
begin real Reone , Iraone , Reother , Imother ; 

i : =0 ; Reone : =one ; Reother : =other ; 

i : si ; I mono : =one ; I mo the r : mother ; 

res : =Reone*Imother+Imone*Reother ; i : =0 ; 

res : =Reone*Reother~Imone* Imother 

end ; 



procedure cd( res, one, other) ; real res, one, other ; 
begin real Reone , Imone , Reother , Imother , denom ; 

i:=0; Reone :=one; Reother : =other ; 

i : =1 ; Imone : =one ; Imother : =other ; 

denom: =Reother*Reother+Imother*Imother ; 

res : =( Imone*Reother -Reone* Imother ) /denom ; 

i : =0 ; res : =( Reone*Reother+ Imone* Imother ) /denom 

end ; 



real procedure real (variable) ; real variable; 
real :=(if i=0 then variable else 0.0); 
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Jggj- £££££^ Wp imaginary ( variable) ; roai variable; 
imaginary ;=( if i=0 thon .0 eigo variable) ; 



real p£CCj3dure cxconj(it) ; real it; 

cxconj: = (if i=0 thon i"t~£lso -it) ; 



roal proc<3dure mod(it) ; real it; 
bogin real Roit ,lmit ; 

i:=0; Reit:=it; i:=1 ; Imit:=it; 
Bod;=sqrt(Reit*Reit+Imit'"Imit) 

end; 



real p rocedu re arg(it) ; real it; 
begin Teal Reit ,Ioit ; 

i:=0; Reit:=it; i:=1 ; Imits=rit; 
aig:=(if Reit>0„0 the n arctan (irait/Reit) els 
if abs(lnit)>0 o then 
sign(lmit)*1 ,57079633 - 
aretan(Reit/Imi t) ) 
else 3 c 1 41 59265 
e end; p* 



ISSUE TWO 
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procedure polar form(res,r ,theta) ; real res,r,theta; 
begin real rl , thetal ; 

rl7=r; thetal :=theta; i:=0; res :=r1*cos< thetal) ; 
; res : =r1 *sin( thetal ) 

end ; 



procedure imult(res,it) ; real res, it ; 
begin real aux; 

i:sO; aux:=it; i:=l; res:=aux; aux:=it; 

1 : =0 ; res:=-aux 

end ; 



procedure compsq(res, it) ; real res, it ; 
begin real Reit, Imit; 

i:=0; Reit:=it; i:=1; Imit:=it; 

res : =2. 0*Reit*Imit ; 

i : =0 ; res : =Reit*Rei t-Imi t* Imi t 

end ; 



procedure comprecip(res,it) ; real res, it; 
begin real Reit t Imit,denom; 

i:=0; Reit:=it; i:=1 ; Imit:=it; 

denom: =Reit *Reit+Imit*Imi t ; 

res : =-Imit/denom; i : =0 ; res : =Reit/denom 

end; 



t 
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procedure cxsqr t<res ,it) ; real res, it; 

poler f crm(res ,sqrt (mod< i t) ) ,0 e 5*arg(it) ) ; 



procedure oompln(res,it) ; real res, it; 
begin real aux; 

aux: = ln(mod(it)) ; i:=0; res— aux; 
aux:=arg(it) ; i:=1 ; res:=aux 

c nd ; 



P ^^^E££ compexi(res,it) ; real res.it; 
b^gin real auxl ,aux2; 

i:=C ; auxl :=exp(it) ; i := 1 ; aux2;=it; 

res :=aux1 *sin(aux2) ; 

i:=C ; res:=aux1 *cos(aux2) 

end ; 
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procedure hyp<sinh,cosh,y) ; value y; real sinh,cosh,y ; 
begin real yl ; 

y1 :=exp(y) ; cosin=0. 5*(y1+1 ,0/yl ) ; 

if abs(y)>no then sinh:=G.5*(y1-1 .O/yl) 

else begin integer r; real br f brplusl , brplus2 ; 

array CWC[0:5]; 

CTC[o]: =1.1 3031 821 ; 

CWC[l]:=4. 43368 498; 

CWC[2]:=5. 42926 312; 

CWCt3]:=3.19843 646j -6; 

CV/c[4]: =1.10367 7 10 -8; 

CWcr5]:=2.498 10 -11; 

brplusl : =brplus2 : =0. ; 

y 1 : =2 . 0*(2 . 0*y*y~l . 0) ; 

for r:=5 step -1 until do 



begin br:=y1*brplusl-brplus2+C 1 iyc£r] ; 

if r£0 then begin brplus2:=brplusl ; 

brplusl :=br 



end 



end ; 

sinh:=y*(br -brplusl ) 



end 



end ; 




procedure compsin(res,it) ; real res, it; 
begin real Reit , Imi t , sinhlmit , coshlmit ; 

i:=0; Reit :=it; i:=l; lmit:=it; 

hyp( sinhlmit , coshlmit , Imit ) ; 

res : =cos< Reit ) * s inhlmit ; i : =0 ; 

res : =sin(Reit)*coshImit 

end ; 



procedure compcos(res,it) ; real res, it; 
begin real Reit, Imit , sinhlmit, coshlmit ; 

i:=0; Reit:=it ; i:=1; Imit:=it; 

hyp( sinhlmit, coshlmit, Imit) ; 

res : =-sin(Rei t ) *sinhlmi t ; i ; =0 ; 

res:=cos(Reit)*coshImit 

end ; 
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tegin real Rei t , Imi t , sinhlmi t , ooshlmi t , sinRei t , cosRei t 



i:=0; Reit:=it; i:=1; Imit:=it; 

hyp(sinImit,coshImit f Imit) ; 
sinReit;=sin(Reit) ; cosReit :=cos(Reit) ; 
auxl [ 1 ] : =cosReit*sinhImit j 
aux2[l ] :=-sinReit*sinhImit ; 

aux itol Z =sinRei t*coshImit ; aux2[0] 
aux2 [ ] : =cosRei t ♦cosh Imit; 
c d ( r es , aux 1 [ i ] , aux2 [ i ] ) 




array auxl; aux2[ :1]; 



end ; 




procedure cxarcsin(res, it) ; real res, it; 
begin real x,y,d1 ,d2,d3,d4; 

i:=0; x:=it; i:=1 ; y:=-it; 

d3:=x*x; d4:=y*y; di:=d3+d4; d2:=d3-d4; 

d3: =di *d1-2. 0*d2 ; d4 : =sqrt (d3+1 . Q) ; 

res : =sign(y) *0 % 5* 

ln( d1 +d4+sqrt ( d3+d 1 * ( dl +2 # 0*d4 ) ) ) ; 

i : =0 ; res : =arct an(x*sqrt (2.0/(1. Q~d2+d4) ) ) 

end; 



procedure cxarccos(res,it) ; real res, it; 
begin array aux£o:i]; 

cxarcsin(aux[i ] , it) ; 

eq(res , real( 1 # 57079 633) -auxti]) 

end ; 



procedure cxarctan(res,it) ; real res, it; 
begin array auxl ,aux2[0: l] ; 

eq(auxl[i],it); compsq( aux2[ i ] , auxl [i]); 

cxsqrt(aux2[i),real(l.0)+aux2[i]) ; 

cd(aux2[i],auxl[i],aux2[i]); 

cxarcsi n( r es , aux2 [ i ] ) 

end ; 
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procedure onehochot her < res , one , other ) ; real res, one, other 
begin array auxl[c:l3; 

compln( auxl C i ] , one) ; em<aux1 [i] , other, auxl [i]) ; 

compexp ( res , auxl til) 

end ; 



boolean procedure even(integer) ; integer integer; 
even: integer =( integer div 2)*2) ; 
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procedure LOGC(a,b,c,d,FAIL) ; value a, b; real a,b,c,d; label FAIL; 
if a=0 and brO then goto FAIL else 
begin real e f f ; 

e:=0.5*a; f:=0.5*b; 

if abs(e) < 3.5 and abs(f) < 3.5 then 
begin c:=abs(2*a)+abs(2*b) ; 
d : =8*a/c*a+8*b/c*b ; 
c : =0 # 5*<ln(c )+ ln(d) ) -1 . 03972077 

end 

else begin c:=abs(c.5*e)+abs<o.5*f ) ; 

d: =0. 5*e/c*e+0. 5*f /c*f ; 

c x =0. 5*<ln(c)+ln(d> ) -f 1 . 03972077 

end ; 

d:=if a $ and abs(e) > abs(f) then arctan(b/a)+(if sign(a) £ -1 then : else 

if sign(b) i -1 then 3.14159265 else -3.14159265) else -arctan(a/b)+1 . 57079633 * sign(b) 

end LOOC; 



SUM SERIES OF CHEBYSHEV POLYNOMIALS 

real procedure Chebsum^PjU^s) ; 
value x,ra,s; real x; integer m,s; array P; 
comment HUGO LIBRARY PROCEDURE FP01 
AUTHOR: W. S. WARM; 

Eva lua te s the s urns : - 
m 

E P[i]T.(x) if s « 1 

m 

£ P[i]T 2 i(x) if s = 2 
1*0 

m 

I P[i]T 2 i+i(x) if s = 3 
i*0 

where T t (x) is the ith Chebyshev polynomial defined by 
Ti(x) = cos(iarcos(x)) , 

The evaluation uses Chenshaw 1 s method involving the 
recurrence relation for Chebyshev polynomials. 

see Chenshaw, C.W. : A note on the summation of 
Chebyshev Series M.T.A.C. 9 188-120 (1955). 



real procedure Chebsum (x,P,m,s) ; 
value x , m , s ; real x; integer m,s ; array P ; 
comment This procedure evaluates the sum 
£# 

> P T (x) where T (x)= cos(i*arcos(x)) and the prime denotes 
i=0 i i i 
halving of the first term, when s=1 . 

m, m 

For s=2 and s=3,the sums > P T (x) and > P T (x) are evaluated 

i=0 i 2i i=0 i 2i+1 

respectively. ; 

beg in real CO, CI ,C2 f y,mult; 

integer i; 

switch SI := LI ,L2,L6; 
switch S2:= L3,L7,L4,L5; 
y:=x; 

goto SI [s] ; 

L2:L6: y:= 2*x*x - 1 .0; 

LI : CI :=C0:=C2:=0 .0; 

mult s= 2 # 0*y; 

tor i:=m £teg -1 until do 

begin C2:= CI ; CI := CO; 

C0:= mult*C1 -C2 + P[i] 

end; 

goto S2[s]; 

L3:L7; Chebsum := # 5 * (CO - C2) ; goto L5; 

14: Chebsum := x* (CO -CI); — — 

L5: end of Chebsum; 



FP01 

J 



LOGARITHM OF FACTORIAL (n) 



real procedure &ogfac(a) ; Integer a; 
comment HUCC LIBRARY PROCEDURE FSOl : 
AUTHOR COMPUTER BULLETIN : 

Evaluates £n(f actorial(a)) . For values of a 4 7 
the procedure computes factorial (a) and uses the 
standard function in to evaluate £n(f actorial(a) ) . 
For values of a > 7 the procedure uses the 
approximation 

£n(f actorial(a)) = &n(gamma(a+l) ) 
* in(/2u) + (a+0.5)*£n(a) - a + Bl/(2a) - B2/(12a) 
where Bl(=l/6) and B2(«l/30) are the first two 
Bernoulli numbers. 



F301 



real procedure logfac(a); integer a; 
comment evaluates ln(f actorial(a)) ; 
begin real aa,k,b; 
aa:=a; 

if aa>7. 5 then begin k:=1 O/aa; 

logfac:=(aa+0 # 5) *ln< aa)+k*<0.0833333333-k*k/360.0) -aa+O. 9189385333 
end 

else begin aa:=aa+0.5; b:=1.0; 
for k:=2 # step 1.0 until aa do b:=b*k; 
logfac:=ln(b) 
end 

end logfac; 
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FOURTH ORDER RUNGE-KUTTA 



procedure RKFOUR(x,y,f ,h,n) ; value h,n; integer n; 
real x,h; array y,f ; 

comment HUCC LIBRARY PROCEDURE GL01: 
AUTHOR: J. BOOTHROYD : 



A procedure for solving n first order linear differential 
equations of the form 

- fi(x,y lf y 2 y n ) 



dyi 
dx 



j[± - f 2 (x,yi,y 2 y n ) 

dx 

^n = f n (x,yi,y 2 Y n ) 

dx 

for known starting values (y^) , (y 2 ) $ (y rt ) 

n 

The correct use of RKFOUR requires the following preliminary 
operations ♦ 

1. Array y must be initialised so that, for y[l:n], 

y[i] = (yJ i « 1,2, n 

1 

2. The right-hand sides of the equations must be specified 
by the declaration of a procedure function such that the 
operation of function places in array f [l:n] the 
appropriate f [i] , i ■ 1,2,. • ..n. 

Examples : A dyj _ 

dx" 

f!i - -k 2 yj 

dx 

procedure function(x,y,f) ; real x; array y ,f ; 

begin x:=x; f[l]:-y[2]; f [2] := -k*k*y [1] end; 

B dy T , 

J 1 = y 2 + x 

dx 

s y 

dx 

fll a -Y2 Y3 - yi x2 
dx 



-2- 
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cont. 



procedure function (x,y ,f ) ; real x; array y,f ; 
begin f[l]:- y[2]+x; f[2]:= y[3]; 

f[3]:= -y[2]*y[3]-y[l]*x*x end; 
A single call of RKFOUR replaces the values y [i] ,i=l,2, . . .n 
appropriate to some x by the corresponding values at x+h. 
The user must, therefore, arrange to call RKFOUR successively 
the appropriate number of times consistent with h, the 
interval of interest and the values of the argument at which 
values of the solution are required. 

Suppose values of y [1] are required at intervals of 0.2 
from x«0 to x=5 inclusive and the computation step interval 
h is 0.05. 

x:= 0; 

for i 1 step 1 until 26 do 
begin print x,y [1] ; 

for j := 1 step 1 until 4 do RKFOUR ( ) 



procedure RKFQUR(x,y,f ,h,n) ; value h,n; integer n; real x,h; array y,f; 

begin real hby2,hby6; integer i ; array ybar,p[l:n] ; 

procedure step(a,b,kl ,k2) ; value k1,k2; real k1,k2; array a,b 

for i:=l step 1 until n do 

begin p[i] :=p[i]+k1*f [i] ; a[i]:=y[i]+k2*b[i] end step; 

hby 2 : =h/2 # ; hby S : =h/6 . ; 

for i:=1 step 1 until n do p[i] :=0.0; 

function(x,y,f); step(ybar,f ,1 # 0,hby2) ; x:=x+hby2; 

f unc t ion(x, ybar , f ) ; step(ybar, f , 2. 0, hby2) ; 

function(x,ybar,f) ; step(ybar,f ,2.0,h) ; x:=x+hby2; 

f unct ion(x, ybar , f ) ; step(y,p,l .0,hbyS) 
end RKFOUR; 



r- c 
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MATRIX INVERSION 

procedure mxinvert(a,n,eps, singular) ; value n,eps ; 

array a; integer n; real eps ; label singular; 
comment HUCC LIBRARY PROCEDURE LE01: 
AUTHOR: J, BOOTHROYD 

Inverts a matrix in its own space using the 

Gauss-Jordan method with complete matrix pivoting. 

I.e., at each stage the pivot has the largest absolute value 

of any element in the remaining matrix. The 

coordinates of the successive matrix pivots used 

at each stage of the reduction are recorded in the 

successive element positions of the row and column 

index vectors r and c. These are later called upon 

by the procedure mxperm which rearranges the rows 

and columns of the matrix. If the matrix is singular 

the procedure exits to an appropriate label in the main 

program. 

This procedure uses procedure mxperm, and for correct 
operation either NC02 or NC04 must previously have 
been declared* 



] LS01 

procedure mxinvert (a, n,eps, singular) ; value n ? eps; array a; i nteger n; real eps; label singular; 
begin integer i, j,k,pivi,pivj,p,ri,ci,rk,c j,iless1 ; real pivot; i nteger array r,c[1 :n] ; 
~ comment sot row and column index vectors; 
fOT~i7=1 step 1 until n do r[i] :~o[i] :=i ; 

comment find initial pivot; pivi:-pivj:=1 ; pivot :=a[l . 1 ] ; 
for i:=1 step 1 until n do for j;=1 Bteg 1 until n £o 

i f 'abs < a [ i , j 1 )>abe (pivot ) then begin pivi:=d; pivjs=j; pivot s=a[i,j3 end; 
comment start reduction; 
for i : =1 step 1 until n do 

begin ri : =rTpivi] ; ripivi] ;=r[i] ; r[i]:=s»i; cis=cEpivjl; c[piv j] :=c[i] ; c[i]:=ci; ilessl :=i-1 ; 
if eps > abs(a[ri,ci]) then 

begin print punoh<3) ,sameline, digits (3) ,££1?MATRIX SINGULAR? 9 C21?i=? , i p ££1?FIV0TS FOLLOIV?; 
±cr i:=1 sten 1 u ntil n do 

print punch( 3) 9 sameiine , digi t s ( 3) * ££1?? , r [ i ] f ££s4?? , c [ i ] ; 

goto singular 

end; 

f or j:=1 step i u ntil ilessi r i+1 steg 1 unt^il n do 

begin c j:=c[ j] ; a[ri,cj] :=a[ri,c j]/pivot end; 
a[ri ,cij:=1. O/pivot ; pivot : =0 ; 
for k:=1 sto p 1 until ilessi, i+1 step 1 unti l n do 
begin rk:=r[k] ; 

for j s=1 step 1 until iles3l ,i+l step 1 until n do 

begin cj:=c[ j] ; a[rk,c j] :=a[rk,c j]-a[ri,c j]*a[rk,ci] ; 

if k>i and j>i and abs(a[rk,c j]) >abs(pivot) then 

begin pivi:=k; pivj:=j; pivot :=a[rk,cj] end conditional 

end jloop; 
a[rk,ci] :=-a[ri,ci]*a[rk,ci] 
end kloop 

end iloop and reduction; ^ ^ J^- 

comment rearrange rows; mxperm(a[ j,p],a[k,p] , j,k,r,c,n ; p) ; 

comment rearrange columns; mxperm(a[p, j],a[p,k], j,k,c,r,n,p) 
end mxinvert; 
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SOLVE LINEAR EQUATIONS - ONE R»H. SIDE 

procedure SOLVEQ(a,n); value n; integer n; real array a; 
comment HUCC LIBRARY PROCEDURE LE02 ; 
AUTHOR: J, BOOTHROYD 

A procedure which solves Ax = b where A is [l:n,l:n] 
and b is [l:n] . The elements of A and b must 
occupy the elements of a[l:n f l:n+l] as follows :- 
a[i,j] = A[i,j] i « 1,2,,.. ,n, j « l,2,.,,.n. 
a[i f n+l] » b[i] i - l,2 >t ».,n f 

i.e, the R.H.S. vector b occupies column n+1 of the 
augmented matrix a. On exit from the procedure the 
solutions x[i] occupy a[i,n+l] i = 1,2,.. t n. 
The method is Gauss-Jordan reduction to diagonal form 
with partial pivoting. No row exchanges are performed 
(the matrix is reduced in-situ) , but a final permutation 
of a[i,n+l] is performed to reorder the solution vector* 
The procedure does not include a test for singularity or 
ill-condition, and should not be used if the equations 
may be suspected of either attribute. 



j LE02 



procedure SQLVEQ(a,n) ; value n; integer n; real array a; 

begin integer i, j,k,piv,ri t rk,n1 ; real pivot, arki; integer array r[l:n]; switch S:=L; 
for i:=1 step 1 until n do r[i]:=i; ni:=n+1 ; 
for i:=l step 1 until n do 

feegin piv : =i ; ri : =r [i ] ; pivot : =a[ri ,i] ; 

. for k:=i+t step 1 until n do 7 
if abs(a[r[k],i])>abs(pivot) then [ u t 

begin piv:=k; pivot : =a[r[k] ,i] end; ) 
ri:=r[piv] ; r[piv]:=r[i] ; r[i]:=ri ; 
for j:=i+1 step t until nl do a[ri , j ] : =a[ri , j]/pivot ; 
for k:=1 step 1 until 1-1*1+1 step 1 until n do 
begin rk:=r[k]; arki :=a[rk,i] ; 

for j:=i+l step i until nl do a[rk, j]:=a[rk, j]-arki*a[ri, j] 
end for k; 
end reduction; 
for i:=n step -1 until 2 do 

begin k:=r[i]; 
L: if k#i then 

begin if K>i then begin k:=:r[k]; goto L end ; 

pivot :=a[i,nl]; a£i,m] :=a[k, nl ] ; a[k,m]:=pivot 

end 

end 

end SCLVEQ; 



j 
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MATRIX DECOMPOSITION Ax = b 

procedure AT0LU1 (a,r,n,eps , singular) ; value n,eps; 

integer n; real array a; Integer array r; label singular; 
comment HUCC LIBRARY PROCEDURE LE03 

AUTHOR: J. BOOTHROYD 

Performs an in situ decomposition of a matrix A into quasi 
triangular matrices L, a lower triangle, and U a strictly 
upper triangle, by Gaussian elemination with partial 
pivoting^ The successive pivots of the reduction constitute 
the quasi diagonal elements of L, These are chosen to be 
the elements of maximum modulus in the leading column of 
successive reduced matrices. The pivotal row subscripts 
are recorded in the non-local vector r[l:n] , the ith 
pivo tal row subscript occupying r[i] t If the matrix is 
singular (eps is the appropriate tolerance) the procedure 
records the ith stage of the reduction at which singularity is 
detected by changing the sign of r[i] before it exists 
to the main program via the label parameter singular. 



i LE03 j 

! t 

1 — --, i , -I , , ™J 

procedure AT0LU1 ( a, r,n,eps f singular) ; value n,eps ; integer n; rfeal eps; 
real array a; integer array r; label singular; 

begin integer i, j,k,piv, ri,rk; real pivot, arki; 
for i:=1 step 1 until n do r[i]:=i ; 
for i:=1 step 1 until n do 

begin pivs=i; ri:=r[i] ; pivot :=a[ri,i] ; 
for k:=i+1 step 1 until n do 

if abs(a[r[k],i]»abs(pivot) then begin piv;=k; pivot:=a[r[k],i] end ; 
if abs(pivot)<eps then 

begin print punch( 3) , ££1?MATRIX IS SINGULAR?; rlpiv] :=-r[piv] ; 
goto singular 

end ; 

ri:=r[piv] ; r[piv] :=r[i] ; r[i]:=ri; 

for j:=i+1 step 1 until n do atri, j]:=a£ri f j]/pivot; 

for k:=i+1 step 1 until n do 

begin rk:=r[k] ; arki :=:a[rk f i]; 

for j:=i+1 step 1 until n do a[rk, j]:=a[rk, j]-arki*a[ri, j] 

end for k 

end for i 

end ATGLU1 ; 
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FORWARD AND BACK SOLUTION Ax = b 

procedure LUlSOL(a,b ,r,n) ; value n; real array a,b; 

integer array r; integer n; 
comment HUCC LIBRARY PROCEDURE LE04 

AUTHOR: J. BOOTHROYD 

Solves the equation Ax~b by the successive operations 
Lf«b and Ux=f where L and U are respectively a lower and 
strictly upper triangle such that L(UH-I)=A. The elements 
of L and U share the array a in accordance with the 
pivoting strategy of procedure AT01U1, 

i.e. L[i, j] occupy a[r[i] , j] j < i i e l,2 ff ,...,n 

and U[i, j] occupy a[r[i],j] ] > i i= 1,2, ,n, 

r [l:n] is a non-local permutation vector in which is stored 
the a array row subscripts of successive rows of L, The 
procedure uses procedure sigma (MS 01) . 

Notes : For a matrix A[l:n,l:n] , vector b[lsn], the procedure 
calls AT0LU1 (A,r ,n, eps ,bad) LE03 
LUlS0L(A,b,r,n) LE04 
PERMB(b,r,n) NC01 
Solve Ax=b for one right hand side* 
For large sets of equations with many right hand 
sides use 

ATOLU 1 ( A , r , ep s t b ad ) LE03 
— Read b 

LUlS0L(A,b,r,n) LE04 
PERMB(b,r,n) NC01 
— Print b 

Alternatively the output procedure ZP01 may be used 
in place of PERMB, Print b. 



procedure LU1 SOL(a,b,r,n) ; value n; real array a,b; integer array r; integer 
comment solves the equation Ax=b by the successive operations Lf=b and Ux=f 
where L and U are respectively a lower and strictly upper triangle such that 
L(U+I)=A. The elements of L and U share the array a in accordance with the 
pivoting strategy of procedure AT0LU1 # 

i.e # L[i, j] occupy afr[i],j] j<i i= 1,2, ,n 

and U[i, j] occupy a[r[i], j] j>i i= 1,2, #### .,n. 
r[l :n] is a non-local permutation vector in which is stored the a array row 
subscripts of successive rows of L. The procedure uses procedure sigma; 
begin integer i, j,ri,ilessl; 
for i : =1 step 1 until n do 

begin ri:=r£i] ; ilessl :=i-i ; 

bfri]:=(btri]-sigma(a[ri,j]*b[r[j]],j,l,ilessl))/a[ri,i} 
end; 

for i:=n-l step -1 until 1 do 
begin ri :=r[i] ; 

bCri]:=b[ri]-sigma(a[ri,j]*b[r[j]],j f i+i,n) 

end 
end LU1 SQL; 



r — ' 
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MATRIC DECOMPOSITION BY GROUT'S METHOD 

procedure crout(a,r ,n, eps , singular) ; value n,eps ; integer n; 
real eps; real array a; integer array r; label singular; 
comment HUCC LIBRARY PROCEDURE LE05 : 
AUTHOR: J. BOOTHROYD : 

Performs an in situ decomposition of a matrix A into quasi 
triangular matrices L, a lower triangle and U, a strictly 
upper triangle, by Croufs method with partial pivoting. 
The successive pivots of the reduction constitute the 
quasi diagonal elements of L. These are chosen to be the 
elements of maximum modulus in the leading column of 
successive reduced matrices . The pivotal row subscripts 
are recorded in the non-local vector r[l:n], the ith 
pivotal row subscript occupying r[i] • If the matrix is 
singular (eps is the appropriate tolerance) the procedure 
records the ith stage of the reduction at which singularity 
is detected by changing the sign of r[i] before it exits to 
the main program via the label parameter singular. This 
procedure uses the procedure sigma. 

The results of this procedure are identical with those of 
LE03. LE05 is not a direct replacement for LE03 since its 
identifier is different and LE05 uses MS01. 



1 



grocedure crout(a,r,n,eps,singuiai-); value n,eps; integer n; rva* eps; 
real array a; integer array r; label singular; 

begin integer i, j,k,piv,ri,rk,iless1 ; real pivot, arki; 
for i:=1 step 1 until n do r[i]:=i; 
for i:=1 step 1 until n do 

begin pivot:=0*0; piv:=i; iless1:= i-1; 
for k:=i step 1 until n do 

begin rk:=r[k]; 

arki:=a[rk,i]:=a[rk,i]-sigma(a[rk t j)*a[r[j],i],j,1,iless1); 

if abs(arki»abs(pivot) then begin piv:=k; pivot:=arki end; 

end ; 

if abs(pivot)<eps then begin print punch(3) ,££1?MATOIX IS SINGULAR?; 

r[piv]:=-r[piv]; goto singular 

end ; 

ri:=r[piv]; r[piv]:=r[i] ; r[i]:=ri; 

for j:=i+1 step 1 until n do a[ri, j]:=<a[ri, j]-sigma<atri,k]*a[r[k] f jljk^ilessDVpivot 
end for i 

end crout; 



TRIANGULAR MATRIX PRODUCT LU=A 



procedure LUlTOA(a,r ,n) ; value n; integer n; real array a; 
integer array r ; 

comment HUCC LIBRARY PROCEDURE LE06: 
AUTHOR: J. BOOTHROYD : 

Performs an in situ matrix multiplication of quasi- triangular 
matrices L, a lower triangle, and U, a strictly upper triangle 
to form a matrix A given, in effect, by A=L(U+1) • The 
elements of L and U are distributed within the array a 
in accordance with the pivoting strategy of procedure AT0LU1, 
i,e. elements L[i, j] occupy positions a[r[i],j] j<i, 1=1,2, • . ,n 
elements U[i , j ] occupy positions a[r[i],j] j>i, 1=1,2, ,.,n. 
r [l:n] is a non-local permutation vector in which is recorded 
the a array row subscripts of successive rows of L. 
Note: LE06 is the inverse of LEO 3, For a matrix A 
the operations 

ATOLU l(A,r, n,eps,s ing ular) 
LUlTOA(A,r,n) 
will leave A unchanged except for corruptions to 
the elements of A caused by computational errors. 



procedure LU1TDA(a,r,n) ; value n; integer n; real array a; integer array r; 
comment performs an in situ matrix multiplication of quasi -triangular matrices L, 
a lower triangle, and U, a strictly upper triangle to form a matrix A given, in 
effect, toy A=L(U+I) . The elements of L and U are distributed v/ithin the array 
a in accordance with the pivoting strategy of procedure ATDLU1, i # e # 
elements L[i, j} occupy positions a[r[il, j] j<i f i=1 ,2, . , . . # . # . , n 

elements U[i, j] occupy positions a[r[i],j] j>i, i=1,2, # .. , n« 

r[l:n] is a non-local permutation vector in which is recorded the a array row 
subscripts of successive rows of L; 
begin integer i» j,k,ii f ri; real sum; 
for i : =n step -1 until 1 do 
begin ii:=i+1; ri:=r[i]; 

for j:=n step -1 until ii do 
begin sum: =a[ri f j]*aCri f i ] ; 

for k: = i-1 step -1 until 1 do sum:=sum+a[ri,k]*a[r[k], j] ; 
a[ri, j]:=sum 
end; 

for j:= i steg -1 until 1 do 
begin sum:=:a[ri f j] ; 

for k:= j-i step -1 until 1 do sum: =sum+a[ri ,k]*a£r[k] , j] ; 
atri, j]:=sum 

end 

end 
end LU1TOA; 
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INVERT QUASI LOWER TRIANGULAR MATRIX 

procedure INVL(a,r,n) ; value n; integer n; real array a; 
integer array r; 

comment HUCC LIBRARY PROCEDURE LE07: 
AUTHOR: J. BOOTKROYD : 

Inverts a quasi triangular matrix L in its own space. 
The elements of the lower triangle L are distributed 
within the array a in accordance with the pivoting 
strategy of procedure AT0LU1, i.e. 

elements L[i,j] occupy positions a[r[i] , j] j<i ,i=»l,2, . , ,n 
r [l:n] is a non-local permutation vector in which is 
recorded the a array row subscripts of successive rows of L. 
The procedure uses procedure sigma (MS01). 
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INVERT QUASI UPPER TRIANGULAR MATRIX 

procedure INVUl(a,r ,n) ; value n; integer n; real array a; 
integer array r ; 

comment HUCC LIBRARY PROCEDURE LE08 : 
AUTHOR: J. BOOTHROYD : 

Inverts a quasi-triangular matrix in its own space. 
The elements of U, a strictly upper triangle, are 
distributed within the array a in accordance with 
the pivoting strategy of procedure AT0LU1, i.e. 

elements U[i,j] occupy positions a[r[i] , j] j>i i=l,2,..,n. 

r [l;n] is a non-local permutation vector in which is 
recorded the a array row subscripts of successive rows 
of U. The procedure uses procedure sigma (MS 01) . 



LE07 
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procedure INVL(a,r «n) : value n; integer n; real array a : integ er ar ra y r ; 

be S i , n integer i ,j ,k .ri .ilessl : real ariiy~~ — =e= 
Jog step 1 until n do" 

begin ri:=r[i]; ilessl :=i-1 ; 
arii:=a[ri,i] ; 

cS2£ J^ 1 1 until ilessl do a[ri,j]:=:-sigma(a[ri,k]*aCr[k],j],k,j f iloss1)/arii 
a[ri,i]:=1 # 0/arii 



end 



end INVL; 
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procedure INVU1 (a,r ,n) ; value n; integer n; real array a; integer array r- 
SSSiSi integer i.i.k.jplual , -jlessl T r>i ; ~ — ^ > 

Jog i:=n-1 step -1 until 1 do 

begin ri:=r[i] ; iplusl :=i+l ; 

|Q£ j:=n step -1 until iplusl do 
begin jlessl :=j-1 • 

a[ri,j]:=-(a[ri,j]+sigma(a[ri,k]*a[r[k3,j],k, iplusl, jlessl)) 

end 

end 

end INVU1 ; 
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TRIANGULAR MATRIX PRODUCT UL=A 

procedure UlLTOA(a,r ,n) ; value n; integer n; integer array r; 
array a; 

comment HUCC LIBRARY PROCEDURE LE09 : 
AUTHOR: J. BOOTHROYD : 

Perfjorms an in-situ matrix product (U+I)L«A where U 
and L are respectively a strictly upper triangle and 
a lower triangle occupying the array a in accordance 
with the pivoting strategy of procedure AT0LU1 

i.e. L[i,j] occupy a[r[i],j] j<i i=l,2,....,n 
and U[i,j] occupy a[r[i],j] j>i 1*1,2, .... ,n. 
r[l;n] is a non-local permutation vector in which is 
recorded the a array row subscripts of successive rows 
The procedure uses procedure sigma (MS01) . 



procedure U1LT0A( a, r t n) ; value n; integer n; integer array r ; array a; 
comment performs an in- situ matrix product (U+I)L=:A where U and L are respectively 
a strictly upper triangle and a lower triangle occupying the array a in accordance 
with the pivoting strategy of procedure ATOLU1 

i.e. L[i, j] occupy a[r[i] f j] j<i i=1 ,2, . , , . ,n 

and U[i f j] occupy a[r[i], j] j>i i=1,2, . 99# ,n, 
r[l :n] is a non-local permutation vector in which is recorded the a array row 
subscripts of successive rows of L. The procedure uses procedure sigma; 
begin integer i, j,k,nlessl,ri; 
nlesst :=n-1 ; 

for i:= 7 step 1 until nlessl do 
begin ri: = r[i]; 

for j:= 1 step 1 until i do a[ri f j]:=a[ri f j]+sigma(a[ri,k]*a[rtk3 f j) f k f i+l f n); 
for j:=i+1 step 1 until n do a[ri, sigma(atri f k]*a[r[k] f j],k, j,n) 

end 

end U1LT0A; 



I LEio] to j LE14 



LIBRARY PROCEDURES LE03 to LE09 produce and perform operations 
on quasi triangular matrices L and Ul where L is a lower triangle 
and U a strictly upper triangle. LIBRARY PROCEDURES LEIO to 
LE14 correspond to some, but not all of LE03 to LE09 for matrices 
LI and U respectively a strictly lower and an upper triangle. 



j LEIO j 

procedure ATOLlU(a,r ,n,eps , singular) ; value n,eps; integer n; 
real eps; integer array r; label singular; 
comment HUCC LIBRARY PROCEDURE LEIO: 

AUTHOR: J. BOOTHROYD : 

see note above and LE03. 



j LE11 

procedure LlUSOL(a,b, r,n) ; value n; real array a,b; 

integer array r; integer n; 

comment HUCC LIBRARY PROCEDURE LE11: 

AUTHOR: J. BOOTHROYD : 

see note above and LE04. 

j LE12 j 

<. i 

procedure INVLl(a,r,n) ; value n; real array a; integer array r; 
comment HUCC LIBRARY PROCEDURE LE12: 

AUTHOR: J. BOOTHROYD : 

see note above and LE07. 



LE10 



to 
contt 



[*LE14 



procedure INVU(a,r,n); value n; integer n; real array a; 
integer array r; 

comment HUCC LIBRARY PROCEDURE LE13: 
AUTHOR: J. BOOTHROYD : 

see note overleaf and LE08. 

LE14 



procedure ULlIOA(a«r«n) ; value nj integer n; array a: 
comment HUCC LIBRARY PROCEDURE LE14: 
AUTHOR: J, BOOTH R)YD : 



see note overleaf and LE09. 



procedure AT0EL1U( a, r,n,eps, singular) ; value n f eps; integer n; real eps; 
real array a; integer array r ; label singular; 
begin integer i f j,k, ri f rk,piv; real pivot, arki; 

for i:=1 step 1 until n do r[i]:=d; 
for i:=l step 1 until n do 

begin piv:=i; pivot :=a[r[i] f i] ; 

for k:=i+l step 1 until n do 

if abs ( a[r[k] , i ] )>abs(pivot ) then begin piv:=k; pivot :=a[r[k] ,i] end; 
if abs (pivot )<eps then 

begin print punch(3),££l?MATRIX IS SINGULAR? ; r[piv] :=-rtpiv] ; 
goto singular 

end ; 

ri : xrtpiv] ; r [piv] : =rf i ] ; r [i ] : =ri ; 
for k:=i+i step 1 until n do 

tegln rk:=r[k] ; arki:=a[rk,i]:=afrk,i]/pivot; 

for j:=i+i steg 1 until n do a[rk f j]:=a[rk, jl-arki*atri, j] 



end 

end 



end ATGL1U; 



i 



I LE11 



procedure L1USOL(a,b,r,n) ; value n; real array a,b; integer array r ; integer n; 
begin integer i, j,ri, ilessl ; 
for i:=2 step 1 until n do 

begin ri:=r[i] ; ilessl :=i~1 ; 

b[ri ] : =b[ ri ] -sigma( a[ri , j]*b[r [ j ] ] , j f 1 ,iless 1 ) ; 
end; 

for i:=n step -1 until 1 do 
begin ri:=r[i]; 

b[ri] : =(b[ri ]«sigma( a[ri , j]*b[r [ j] ] , j, i+1 ,n) )/a[ri , i ] 

end 
end L1US0L; 



I LE12 



procedure INVLl(a,r,n) ; value n; integer n; real array a; integer array r; 
bQ g in integer i , ,j . k . ri . iless 1 . iplusl ; 

for i:=2 step 1 until n do 

begin ri:=r[i] ; ilessl :=i -l • 

for j:=1 step 1 until ilessl do 
begin jplusl:=j+i ; 

atri » j] : =^<a[ri , j]+sigma(a[ri,k]*a[r[k) , j] ,k, jplusl , ilessl ) ) 

end 

end 

end INVL1 ; 



) 



) 



LE13 



procedure INVU(a,r f n) ; value n; integer n; real array a; integer array r; 
begin integer i, j t k,ri,iplus1 ; real arri; 

for i:=n step -1 until 1 do 

begin ri:=r[i] ; iplusl :=i+1 ; arri:=a[ri,i] ; 

for j:=n step -1 until iplusl do atri f j] : =-sigma( a[ri ,k]*a[r [k] , j] ,k, iplusl , j)/arri ; 
a[ri,i] : =1.0/ arri 

end 

endlNVU; 



LEU 



procedure UL1TOA(a,r f n) ; value n; integer n; integer array r; array a; 
begin integer i, j,k,nlessl,iless1 ,ri ; 
nlessl:-n-1 ; 

for i:=1 step 1 until n do 

begin ri:=r[i] ; iless1;=i-1; 

for j:=1 step 1 until ilessl do a[ri, j]:=sigma(a[ri,k]*a[r[k], j] # k,i,n) ; 

for j:=i step 1 until nlessl do a[ri , j] : =a[ri , j]+sigma( a[ri f k]*a[r[k] , j) ,k, j+1 , n) 

end 

end UL1TOA; 



LE15 



SOLVE TRIDIAGONAL LINEAR EQUATIONS 



procedure tridiag(a,b 5 c,d,£o,hi ,eps , f ail) ; value io ,hi ,eps ; 



integer £o,hi; real eps ; label fail 
comment HUCC LIBRARY PROCEDURE LE15 
AUTHOR J. BOOTHROYD 



real array a,b,c,d; 



Solves the equations 
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by the recurrence formulae 



x - d T 
n n 



x ~ d ? - c ' 



r r+1 



di 9 = dx/bx 



d 1 = d -a d - f 
r r r r-1 



b -a c i 
r r r-1 



b -a c - 1 
r r r-1 



r=2,3, ,n. 



The coefficients occupy the elements of arrays a,b,c,d as follows: 
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eps is a singularity tolerance. If l> -a c i' < eps for any r the set 

is presumed singular and exit occurs through the label parameter fail. 

The values of b -a c - f and r are recorded on the control typewriter, 
r r r-1 

The solution is delivered in array d and array c is overwritten in the process; 



procedure tridiag(a,b,c f d, lo f hi, eps,f ail) ; value lo,hi,eps; integer lo,hi; 
real ops; label fail; real array a,b,c,d; 

begin integer i,iless1; real w,ailessl; 

dflol ; =d[lo]/b[lo] ; w;=b[lo] ; 

for i : =lo+1 step 1 until hi do 

begin ilesst :=i-1 ; ailessl :=a[ilessl] • 
cfilessl] : =c[ilessl]/w; 
w:=b[i]-ailessl*ctilessl]; 

if vKeps then begin print punch(3) , scaled( 3) , digits( 3) , 

££1?TRIDIAG MATRIX SINGULAR 
goto fail 

end ; 

dri]:=(dli)-ailesst*d[iless1])/w 

end ; 

for i:=hi-l step -1 until lo do d[i]:=d[i]-c[i]*d[i+l J 
end tridiag; 



s 



I LE15 



-?,w,£Et?i=?,i; 



I 



/ LE 16 I 

DECOMPOSITE SYMMETR IC POSIT IVE DEFINITE MATRIX 

procedur e choleski (a ? n,f ail) ; value n; Integer n; array a; label fr.il; 

comment HUCC LIBRARY PROCEDURE LE 16 (revised) 
Author ; J, Boothroyd. 

Performs the decomposition of a symmetric matrix A into lower 
and upper triangles L,L T such that LLT=A. The elements of the 
lower triangle of A are replaced by the elements of L. The 
strictly upper triangle of A is not used. If the decomposition 
fails, indicating that A is non positive definite the procedure 
exits through the label fail. The procedure uses procedure 
sigma (MS 01); . .. 

LE 17 

INVERT LOWER TRIANGULAR., MATRIX 

procedure linv (a,n) 5 value n; integer n; arr ay a; 
comment HUCC LIBRARY PROCEDURE LE 17 (revised) 
Author : J, Boothroyd 

Replaces the lower triangular portion L of a[l:n,l:n] by the 
elements of L" 1 . The procedure uses MS 01; 

I w ^* 

LE 18 

MATRIX MPLTI PLIGAT1 n T 
— - - - 

procedure mxmult (a ? b ? c,m 5 n,p) ; value m,n,p; 

integer n,n,p ; arr ay a,b,c; 

comment HUCC LIBRARY PROCEDURE LE 18 

Author : J. Boothroyd. 

c[l:m,l:p] 5- a[lm r lsn] * b[l;n ? l:p]; 

This procedure uses MS 01. 



procedure choleski( a,n, fail) ; value i nteger n; array a; label fail ; 
begin integer i , j,k, iless 1 ; real aii; 

for i:=r 1 step 1 until n do 
begin ilessl := 

aii:=a[i,i]-sigma(a[i,k]*a[i,k],k, 1, ilessl) ; 

if aii<j then goto fail; 

aii :=a[i F i] :=sqrt(aii) ; 

for j:= i+1 step 1 until n do 

a[j,i]:= (at j, i] -sigraa(a[i ,k]*a[ j,k] ,k, T ,ilessl ) )/aii 



end 



endcholeski ; 



procedure linv(a,n) • value n; integer n; array a; 
comment inverts the lower triangle of a in situ; 
begin integer i , j ,k, ilessl ; real aii; 
for i:=l step 1 until n do 

begin ilessl:^ aii :=a[i ,i] ; 

for j: = 1 step 1 until ilessl do 

afi, j]:= -sigma(a[i,k]|*a[k, j] f k, j, ilessl )/aii; 
a[i,i]:= 1 # 0/aii 



end 



LE16 



LE17 



end linv; 



procedure raxmult(a, b f c,m, n,p) ; value m,n,p; integer m,n,p; ar ray a,b,c; 
comment c[m pit- a[m n] b[n p] ; 
begin integer i, j,k; 

for i : = 1 step 1 until m do 

for j:= 1 step i until p do c[i , j] : =sigma(a[i ,k]*b[k t j] ,k, T ,n) 



endmxmult ; 



LE18 



LE 19 



DECOMPOSE BANDNATRIX INTO LOWER AND UPPER T TRIANGLES y 
procedure bandlr (a,n,m,eps , escape) ; v alue n,m; 
integer n,m; array a; real eps; label escape; 

comment HUCC LIBRARY PROCEDURE LE 19 
Author : J. Boothroyd. 

Decomposes the n*"* 1 order hand matrix of width m stored in 
a[l:n,l:m] into a lower 'triangle 1 L and a strictly upper 
triangle U. At entry to the procedure the elements of band 
matrix A[l:n,l:n] must occupy the matrix a[l:n,l:m] in 
accordance with the following rules 

A[i f j] occupies a[i,j] l<i< (m-fl) div2 

A[i,j] occupies a[i , j-i+(nrt-l) div2] (irtfl) div2<i<n 

Diagrammatical ly 
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Array A[l:n,l:n] array a[l:n,l:m] 

UNUSED ELEMENTS of matrix a must be set to zero . 
This procedure does not incorporate a pivoting strategy, 
eps is a tolerance which provides a criterion for deciding 
whether any leading submatrix is singular in which event the 
procedure exits to the label escape; 



LE 20 

FORWARD SOLUTION MP BACK SUBSTITUTION FOR BAND EQUATION Ax=b 
procedure bandsol(a,n 3 mb) ; value n,m; 
integer n,m; array a,b ; 

comment HUCC LIBRARY PROCEDURE LE 20 
Author : J. Boothroyd 

Solves the equation Ax=b where A is a bandmatrix by the 
successive operations Lf=b and Ux=f where L and U are 
respectively a lower triangle and a unit upper triangle 
occupying the array a[l:n;l:m] following the use of procedure 
LE 19. The solution vector replaces the vector b[l:n]. 

Procedures LE 19 and LE 20 may together be used to solve 
Ax=b for unlimited right-hand sides by some scheme similar to 

Read Matrix A 

bandlr (A ,n 9 m,eps , escape) LE 19 

^Read vector b 

bandsol (A,n,m,b) LE 20 

Print solution vector b 
I 



LE 19 



procedure bahdlr( a, n,m,eps, escape) ; value n,m; integer n,m; 
begin integer i , j, k t jlim,klim,nless1 ; real pivot; 

jlim:=klim:=(m+1 > div 2; nleasl :=n-l ; 

for i:=1 step 1 until nless 1 do 



real array a; label escape; real eps; 



begin pivot :=a[i, l] ; if abs(pivot )<eps then goto escape; 
for j:=2 step 1 until jlim do a[i , j] : =a[i , j]/pivot ; 
for k:=i+1 step 1 until klim do 
begin pivot : =a[k, f] ; 

for j:=2 step 1 until jlim do a[k, j-1 ] : =a[k, j]-pivot*a[i , j] ; 
' for J- =Jli^ +1 step 1 until m do afk, j-1 ] : =a[k, j] ; 
a[k,m] :=pivot 

end; 

if klimin then klim;=klim+l 

end 



end; 



procedure bandsol(a,n,m,b) ; value n,m; integer n,m; real array a,b; LE 20 

begin integer i , j ,k # lim, jlim; real sum; t 

lim:=(m+1) div 2; jlim:=1; 
for i;=1 step 1 until n do 

begin sum:=:0.0; k:=0; 

fog j:=i~1 £teg -1 until jlim do 

begin snm:=sum+b[ j]*a[i ,m-k] ; k:=k+1 end ; 
b[i] :=(b[i]-Gum)/a[i, i] ; 
if i-jlim+2>lim then jlim:=:jlim+1 
end forward solution; 
jlim:=n; 

for i:=h-l step -1 until 1 do 
begin sum: =0.0; k:=2; 

for j:=i+1 step 1 unti l jlim do 

begin sum:=sum+b[ j]*a[i,k] ; k:=k+l end ; 
b[i] :=b[i]-sum; 

if jlim+2-i>lim then jlim: =jlim-1 
end back substitution 

end bandsol ; 



• LE 21 

SOLVE BAND EQUATIONS Ax=b 

procedure bandmx (a,n,m,b,eps , singular) ; value n,m,eps ; 
real eps ; integer m 3 n; label singular; array a,b ; 

comment HUCC LIBRARY PROCEDURE LE 24 
Author : J# Boothroyd 

Solves the n*"* 1 order linear equation Ax=b where the elements 
of A form, a band matrix of width nw The band elements of 
A[l:nJ.:n] must occupy the array a[l;n,l:m] in accordance 
with the rules 

A[i , j ] occupies a[i , j ] 1< i< (m+1) div2 

A[i,j] occupies a[i ^ i-i+(m+l) div 2] (m+1) div 2<i<n 

Diagrammatically 
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Array A[l:n,l :n] array a[l:n,l:m] 



UNUSED ELEMENTS of matrix a must be set to zero . 
The procedure uses maximum cc lumn pivots and the solution 
vector replaces the elements of vector b[l:n]. eps is a 
singularity tolerance and the procedure exits to the label 
singular if the matrix is singular or excessively ill conditioned; 



LE 21 



procedure ban^xCa>n,m,b,eps ? singular) ; valuo n,m,eps; real eps; integer m,n; ;i|bol singular; ^ a,b; 

feilS integer i , j ,k,raess1 ,liin,piv,ri , rk, iless 1 ; real pivot ,bri , ark 1 , sum; xntege r arrg£ r[l,n] t 

- — — foTi:=1 step 1 until n do r[i]:=i ; nlessl : =n-l ; lim:=:<m+1) djv 2; 

"for i:=T step 1 until nlessl do 

begin piv:=i; pivot : =a[r[i] , 1] ; 

for k:=i+1 step 1 until lim do 

iT"abs<a[r[kl,l])>abs(pivotrthon b3gin piv:=k; pivot :=a[r[h] , 1 ] end; 

if abs(pivot)<eps then begin print punoh(3), ££1?MATRIX SINGULAR?; goto singular end; 
7i:=r[piv]; r[piv]:=r[i] ; r[i]:=ri; 
bri :=b[ri] :=b[ri] /pivot ; 

for j:=2 step 1 until m do a[ri , j] :=a[ri , j]/pivot ; 
for k:-i+i step 1 until lim do 

begin rk:=r['k]; arkl ;=a[rk, l] ; 
b[rk]:=b[rk]-ark1*bri ; 

for j:=2 step 1 until m do a[rk, 3 : =a[rk, j]-arkl *a[ri , j] ; 
"a[rk,m] :=0,0 

end; 

if limin then lim:=lim+1 

end ; 

rit=r[n] ; b[ri] : =b[ri]/a[ri , 1 ] ; 
lim:=2; 

for i:=hless1 step -1 until 1 do 

begin ilessl ; suiii:=O t O; ri:=r[i]; 

for j ; -2 step 1 until lim do sum: =suin+a[ri , j]*b[r[iless 1+ jj ] ; 
b[ri] : =b[ri] -sum ; If limim then lim:=lim+1 

* end; 

for i:=1 step 1 until n do a[i , 1 3 : =b[r[i] ] ; 
for i:=1 step 1 until n do b[i]:~a[i,l] 

end foandmx ; 



LE 22 

REARRANGE PERMUTED TRIANGULAR MATRIX 
procedure MOVE"(a)to: (b) order : (n) pivots : (r) ; value n; 
intege r n ; array a,b ; integer array r ; 
comment HUCC LIBRARY PROCEDURE LE 22 
Author : J, Boothroyd 

LIBRARY PROCEDURE LE 10 decomposes a matrix a into quasi 
triangular interlocking matrices L (a unit lower triangle) 
and 7, an upper triangle which occupy the array a. This 
procedure separates L and \ replacing the elements of U 
formerly in a by zeroes and rearranging D so that it is a 
proper upper triangular matrix in array b • The permuted 
'diagonal' elements of L are set equal to one; 



LE 23 

procedure MOVEL(a) to: (b) order: (n) pivots : (r) ; value n; 
integer n; array a,b; integer array r; 

comment HUCC LIBRARY PROCEDURE LE 23 
Author : J, Boothroyd 

LIBRARY PROCEDURE LE 03 decomposes a matrix a into quasi 
triangular interlocking matrices L, a lower triangle and U, 
a unit upper triangle , which together occupy the array a. 
This procedure separates L and U 9 replacing the elements of 
L formerly in array a by zeroes and rearranging L so that it 
becomes a proper triangular matrix in array b. 
The permuted 'diagonal 1 elements of T T are at equal to one; 



I r— 

LE 22 



procedure MOVElK a) to : < b) order : ( n) pi vot s : ( r ) ; value n; 
integer h; array a,b; integer array r ; 
begin integer i,ri; 

for i:=1 step 1 until n do 
begin ri :=r[i] ; 

b[i,i] :=a[ri f i] ; a[ri , i] : -1 # 0; 
for j:=i+1 step 1 until n do 

begin b[i , j] : =a[ri , j] ; b[ j,i]:=a[ri, j] :=0.0 end 

end 

endM OVE ; 



procedure MDVEL(a) to: (b) order : (n)pivots: (r) ; value n; 
integer h; " array a,b; integer array r ; 
begin integer i,ri ; 

for i :=1 step 1 until n do 
begin ri :=r[i] ; 

b[i,i]:=a[ri,i] ; a[ri,i] :=1 .0; 
for j:=i-1 step -1 until 1 do 

begin b[i, j] :=a[ri, j] ; b[j ,i ] :=a[ri , j] :=0«0 end 

end 

end MOVE; 



LE 23 



LE 24 



DECOMPOSE SYMMETRIC. POSITIVE DEFINITE MATRIX 

procedure SYMDET (a ,n f symde t , f ai 1) ; value n; 
Integer n; real symde t; array a; label fail; 



comment HUCC LIBRARY PROCEDURE LE 24; 
AUTHOR J. BOOTHROYD: < 



Procedure choleski (LE 16) performs the symmetric decomposition 
of a symmetric matrix stored in conventional form. This 
procedure performs an in-situ decomposition of the upper 
triangle of a symmetric, .matrix stored in vector form, 
permitting the solution, in Algol , of symmetric equations 
of order 100 in a 503 with 8K store. The time is approximately 
90 seconds* 

TJie upper triangle of A[l:n,l:n] is stored in a [ 1 : n* (n-f 1) div2 ] 
by columns . 

i.e. The upper triangle of L %A11 A 12 A13 A14 

A2t V A22 A23 A24 
A31 A3V.A33 A34 



is stored as :~ 



A41 A42 A43-.A44 

V 

N 



a All A12 A22 A13 A23 A33 A14 A24 A34 A44 



and element Ai,j is referenced by a[i+i*( j-1) div2 ] . 
To avoid repeated evaluation of subscripts of this form it 
is preferable to set up a column index vector c[l:n] such 
that :=1*(i~l) div2 « permitting references to Ai,j as 
a[i+c[j ] ] . This procedure assumes the exi stance of such a 
global vector which is best Initialised by a recursion relation 
before entry to the procedure using :- 

r.-O; for j :=1 step 1 until n jdo begin c[j] :=r;r:=r+j end 

If the procedure is non-positive-definite exit to the label 
fail occurs . On normal exit from the procedure symdet yields 
the de terminant of the matrix A; 



LE 25 



FORWARD MP BACK SOLUTION Ax=b . SYMMETRIC A 

procedure SYMSOL(a,b ,n) ; value n; integer n; array a,b ; 

comment HUCC LIBRARY PROCEDURE LE 25: 
AUTHOR J. BOOTHROYD: 

Solves the equation Ax=b for symmetric A[l:n,l:n] by the successive 

XT T 
operations Lf=b ,L x=f where A=LL . The triangular matrix L should 

occupy the vector a[l:n*(n-;-l) div 2] using column storage, as it will 

if generated by LE 24, in the description of which procedure find 

details of the method of storage. The procedure assumes the existance 

of a global array c [l:n] initialised before entry to the procedure 

T 

so that c [ j ] :=\i* ( div2, permitting access to an element L i , j by 
reference to a[i+c[ j] ] , On exit froii-: the procedure the elements 
of the solution vector occupy b [1 :n] ; 



( 



pxppedure SYI<IDET(a,n,symdet ? fail) ; value n; integer n; real symdet ; array a; label fail; 
begj-n integer i,j 9 k ? iless1 ,ci 9 ii ,c j ,i j ; 

real det ? aii.aki,aij ; 

dot:=l t 0; 

for is=1 step 1 until n do 

jjog in ilossl :=i-1 ; ci :=c[i] ; ii:=i+ci; aii:=a[ii] ; 
for k:=1 step 1 until ilessl do 

begin aki :=a[k+ci] ; aii :=aii-aki*aki end ; 
If aii<0 *0 then goto fail; 
det:=dot*aii ; 
aii:=a[ii] :=sqrt(aii) ; 
for j:=ri+1 step 1 until n do 

begin cj:=c[ j] ; ij:=i+cj; aij:=a[ij]; 

for k:=1 step 1 until ilessl do ai j :=ai j-a[k+ci]*a[k+c j ] ; 
a[ij]:=aij/aii 

end j 

end i ; 
syiiidet :r:det 
end dYMJET; 



LS 24 



trocedure SYMSC L(a,b,n) ; value n; integer n; real array a,b; 
begin integer i,j, jlessl .cj ; real bi,bj ; 

for j:=1 step 1 until n do 

begin bj:=b[j]; jlessl ;=j-1 ; cj:=e[j]; 

for i:=1 step 1 until jlessl dto bj :=bj-a[i+c j]*b[i] ; 
b[j]:=bj/a[j+cj] 

end ; 

fcr i :=n step -1 until 1 do 
begin Eis=bITjy 

for j: = i+1 ste p 1 until n do bi s=bi-a[i+e[ j]]*b[ j] ; 
b[i]:=bi/a[i+c[i]] 

end 

end SYMSOL; 
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LE 26 

DECOMPOSE BAND MATRIX WITH PIVOTING 

procedure piv2rband(a,n,m s r,eps , singular) ; value n,m,eps; 

n,m; real eps ; array a; label singular; intege r array r; 

comment HUGO LIBRARY PROCEDURE LE 26: 
AUTHOR J. BOOTHROYD: 

Decomposes the n th order bandmatrix of width m stored in columns l:m of 

atl:n»l:mh qdly2 3 into a lower 1 triangle'L and a strictly upper 'triangle 

At entry to the procedure the elements of band matrix A[l:n,l:n] must 

occupy the first m columns of a in accordance with the rules :- 

A[i,j] occupies a[i, j] l<i< (m+1) div 2 

A[i, j ] occupies a[ij j-i+(m-H.) div2] (m+1) div2<i<n 

D i ag r annua ti c a 1 ly : « 

All A12 A13 All A12 A13 

A21 A22 A23 A24 A21 A22 A23 A24 

A31 A32 A33 A34 A35 A31 A32 A33 A34 A35 

A42 A43 A44 A45 A46 A42 A43 A44 A45 A46 

A53 A54 A55 A56 A53 A54 A55 A56 

A64 A65 A66 A64 A65 A66 

Array A[l:n,l:n] Array a[ 1 :n ,1 :m+mdiv2 ] 

jMgSD^LEjffiNTS OF .MATRIX a should be set to zero. 

On exit from the procedure the elements of L occupy column 1 and the 
last m div2 columns of a. The elements of U (emitting the unit 
diagonal elements) occupy the remainder of matrix a« The procedure 
utilises column pivoting without row interchanges with two consequences 
(a) the algorithm is fast and (b) the storage scheme for elements of L 
is complicated. All information relating to the succession of pivotal 
rows is however retained in the index vector r[l:n] for subsequent use 
by LE 27. If the matrix is singular, eps is the tolerance, the procedure 
exits via the label parameter; 
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pivlrband(a,n,m,r , eps , singular) ; value n,m,eps; jjiteger n,m; 
real eps; lg.bel singular; arrag; a; i nteger array r; 

Jbcgin i nteger i, j ,k f p,ri,rk,piv,nloss1 ,lim; real pivot, arkl ; 

fcr i:=1 step 1 until n do r[i]:=i; 
lin:=(nH-1 ) div 2; nlessl :=n-1 ; 
fcjr i:=1 Btgg 1 until nlessl do 

begin piv:=i; ri:=r[i]; pivot :=a[ri ,1 ] ; 
for k:=i+1 s t ep 1 until lim do 

if abs(a[r[k]/l ]>>absTpivot) then begin piv:=k; pivot :=a[r[k] ,1 ] end ; 

if abs(pivot)<eps then begin print puneh(3) ,££1?MATRIX SINGULAR?; goto singular end ; 

if pivfri then b egin ri:=r[piv]; r[piv] :=r[i] ; r[i]:=ri end; 

for jt= 2 s t ep -\ until m do a[ri, j]:=a[ri, j]/pivot; 

p:=m; 

for k:=i+1 s t op 1 until lim do 

begin p:=p+1 ; rk:-r[k]; arkl s=a[ri,p]:=a[rk,1 ]; 

j for j:=2 step 1 until m do a[rk, j-1 ]:=a[rk, j]-ark1 *a[ri, j] ; 
a[rk,m]:=0 # 

end ; 

if liio|=n then lim:=linM-1 

end 

end pivlrband; 



LE 27 

FORWARD AND BACK SOLUTION Ax=b, BANDMATRICES 

procedure pivsolband(a,b ,n,m,r) ; value m,n; integer m,n; 
array a 5 b; integer array r; 

comment HUCC LIBRARY PROCEDURE LE 27: 
AUTHOR J. BOOTHROYD: 



Solves the equation Ax«b where A is a bandmatrix of order n f width m 
by the successive operations Lf=b, Ux=f where the 1 triangular* matrices 
L and U, such that LU^A, occupy the array a[l:n,l:nri-mdiv2] following 
the use of the decomposition procedure LE 26, r[l:n] is a pivotal 
index vector whose elements are generated by LE 26. These are 
subsequently used by LE 27 to effect the correct sequence of operations 
in the forward solution and back substitution; 
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procedurs pivsolband(a,b,n,m,r) ; value zn,n; array a,b; integer array r; integer n,m; 
begin ^integer i , j ,k,ri ,rk, lim; real bri ; 

liia:=(Ei+1) div 2; 

begin ^nt eg or array row[1 :n] ; switch s:=scan; 
for i:=1 step 1 until n do row[i] :=i ; 
for i:=1 stop 1 until n do 

begin ri:=r[i] ; j:=x; 

' if ri£row[i] then begin scan: j:=j+1 ; if row[ j]£ ri then goto scan; 

rov;[ j] :=row[i] ; row[i]:=yi 
^end; ^ 
bri:=b[ri]:=b[ri]/a[ri,1 ] ; 
j:=m; 

I for k:=i+1 stop 1 until lim do 

\ begin ; rk:=row[k] ; b[rk] :=b[rk]-bri*a[ri , j] end; 

if liia4= n then lim:=liEH-1 

end 

end forward solution; 

roal arrajr bi [1 : n] ; 

bi[n]:=b[r[n]]; 

for i:=n-l step -1 until 1 do 

begin ri:=r[i]; bri:=b[ri] ; 
k:=r1 ; 

for j:=i+1 stop 1 until lim do 

begin krrrk+l"; bri:=bri-b[r[ j]]*a[ri,k] end; 
bi[i] :=b[ri] :=bri ; 
if k=m then lim: = liin-1 

end ; 

for i;=1 sto£ 1 until n do b[i]:=bi[i] 
end backsubstitution 

ond v A vcolbanc ; 



LE 28 



DECOMPOSE SYMMETRIC POSITIVE DEFINITE MATRIX 

procedure SYMDET(a,n,symdet,f ail) ; value n; integer n; 
real symdet; array a; label fail; 

comment HUCC LIBRARY PROCEDURE LE 28: 
AUTHOR J. BOOTHROYD: 

Performs the same function as LE 24, using the same method of storing 
the elements of the upper triangle of A[l:n,l:n] as elements of 
a[l:n*(n+l) div2 ] * This procedure does not . however « use the global 
addressing vector c[l:n] necessary with LE 24. It is approximately 
7% faster than LE 24. 
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FORWARD AND BACKWARD SOLUTION OF Ax=b« SYMMETRIC 

procedure SYMSOL (a,b,n); value n; integer n: array a«b: 

comment HUCC LIBRARY PROCEDURE LE 29: 
AUTHOR J. BOOTHROYD: 

Performs the same fraction as LE 25 but avoids 
the need for the addressing vector c[l:n]; 



prcc a duro GYM ET(a , n,symdot , fail) ; value n; integer n; jreal symdet ; arrajr a 
begiL. integer i , j ,ki ,ii ,k1 ,kiloDSl ,i j ,kj ; 

real dot ,aii ; aki ,ai j ; 
d€t: = i .0 ; ii:=1 ; 
for i:~1 step 1 u ntil n do 

begin aii :=a[ii] ; kl :=ii-i+1 ; kilessl :=ii-1 ; 
for ki:=kl step 1 until kilessl do 

begin aki:=a[ki] ; aii :=aii-aki*aki end ; 
if aii<0 then goto fail ; 
det :=det*aii ; 
aii:=a[ii] :=sqrt(aii) ; 
i j :=ii+i ; 

for j:=i+1 step 1 until n do 

begin aij;=a[ij] ; kj:=ij-i+1; 

for ki:=k1 step 1 until kilessl do 

begin ai j :=ai j-afkij^aEkj] ; k j ;=kj+1 end 
a[i j] :=ai j/aii ; i j :=i j+j 

end ; 
ii:=ii+i+1 ; 

end; 
symiet :=rdet 
end SYKDET ; • 



procedu re SYB/I5 fV L<a,b,n) ; value n; integer n; array a,b; 
begin ir^tegor i , j ,k,ii ; real sum; 
ii:=1 ; 

for i:=1 step 1 until n do 
begin k: = ii-1 ; suni:=b[i]; 
toy j:=i-l step -1 until 1 do 

bogin sum:=sum-b[ j]*a[k] ; k:=k-1 end; 

*FfTTT=f3um/a[ii] ; ii :=ii+i+1 
onu forward solution; 
k:-ii-l ; ii:=ii-n-1 ; 
for i:=rn step -1 until 1 do 
beg ill sum:=b[i]; 

for j:=i+1 step 1 until n do 

begin sum:=sum-b[ j]*a[k] ; k:=k+j end; 

bLi] :=sum/a[ii] ; ks=ii-1 ; ii:=ii~i 
end back substitution 
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REDUCE MATRIX TO UPPER HESSENBERG FOBM 

procedure hessenberg(a,n) ; value n; integer n; 

comment HUCC LIBRARY PROCEDURE IJ&30: 
AUTHOR : J. BOOTIIROYD: 

Transforms a real matrix a[l:n,l:n] to upper lies senb erg 
form by a sequence of elementary similarity transformations, 
The transformation is performed using a pivoting strategy 
which involves rot* and column interchanges. The method is 
similar to Gaussian elimination. For a compact CROUT-like 
method which converts a matrix to lower Res senb erg form 
see procedure TRINGLE in Parlett's paper M Languerre f s 
Method Applied to the Matrix Eigenvalue Problem" Math, 
Comp vl8 1964 469 - 485; 



procedure hessenbergC a , n) ; value n; integer n; array a; 
begin integer i, j,k,iless1 ,iplus1 ,piv, nlessl ; real pivot , ink , aki , temp ; 
array m[3:n] ; boolean nonzero; switch s:=pivsearch; 
pivot : =0.0; nlessl :=n~1 ; 
for i:=1 step 1 until nlessl do 
begin iplusl : =i+1 ; 

if pivot =0.0 then goto pivsearch; 
if i4=piv then 

begin for j:=iless1 step 1 until n do begin temp:=a[i, j] ; a[i, j] :=a[piv, j] ; a[piv, j] :=temp end; 
for k:=1 step 1 until n do begin temp:=a[k,i] ; a[k, i] :=a[k,piv] ; a[k,piv] :=temp end 

end ; 

nonzero := false ; 

for k:=iplus1 step 1 until n do 

begin mk : =m[ k] : =a[ k , i 1 ess 1 ] /pivot ; 

if mk^O.O then begin nonzero: = true; a[k,iless1 ] :=0.0 end 

end ; 

if nonzero then 

begin for k:=iplus1 step 1 until n do 
begin mk:=m[k]; 

for j:=i step 1 until n do a[k, j] :=a[k, j]-mk*a[i f j] 

end ; 

piv:=iplus1 ; pivot: =0.0; 
for k:=1 step 1 until n do 
begin aki:=a[k,i] ; 

for j:=iplusi step 1 until n do aki :=aki+m[ j]*a[k, j] ; 

a[k,i] :=aki ; 
if k>i then 

begin if abs(aki)>abs(pivot) then begin piv:=k; pivot :=aki end end 

end 

end 

else 

pivsearch: begin piv:=iplus1 ; pivot : =a [pi v, i] ; 

for k:=i+2 step 1 until n do 
begin aki:=a[k,i] ; 

if abs(aki)>abs(pivot) then begin piv:=k; pivot :=aki end 

end 

end; 

ilessl :=i 

end i 
end hessenberg; 



) 
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TRIDIAGQUALI ZE SYMMETRIC MATRIX GIMS METHOD 
procedure sivens(a,s ,n) ; value n; integer n; array a,s; 

comment HUCC LIBRARY PROCEDURE LE31: 
AUTHOR : J. ROOTHROYD: 

Performs a similarity transformation of a full symmetric 
matrix a[l:n,l:n] to tridiagonal form by GIVENS Method. 
The orthogonal rotation matrices used have the typical 
form: 

10 

Pr = 

c s 
10 

s -c 

-1 t 

so that P = P = P . 

r r r 

The continued product P^,?2 ,P^»,#*of transformation 
matrices is finally available in array s[l:n,l:n] ; 
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procedure givens(a,s,n) ; value n; integer n; array a,s; 
begin integer i,p,q, nlessl ,k; 

real akp , akq, app, apq, aqq, aip, aiq, f ac , cos , ein, cs ; 

for i:=1 step 1 until n do 

begin s[i/ij:=l # 0; 

" for JJ =i +1 _gtep 1 until n do s[i, j] : =s[ j,i] :=0«0 

end ; 

nlessl :=n-1 ; 

for p:=2 step 1 until nlessl do 
begin k:=p-1 ; akp:=a[k,p]; 

for q:=p+1 step 1 until n do 
begin akq:=a[k,q] ; 

" apq:=a[p,q]; app:=atp,p]; aqq:=a[q,q]; 
fac;=sqrt(akp*akp+akq*akq) ; 
if fac*O t O then 

begin cos:=akp/fac; sin:=akq/f ac ; 
for i : =1 step 1 until n do 
begin if i>p then 

begin aip:=a[i,p]; aiq:=a[i,q]; 

a[i,p] :=a[p,i] : =cos*aip+sin*aiq 
a[i,q] :=a[q,i] :=sin*aip-cos*aiq 

end; 

aip:=s[i,p]; aiq:=s[i,q]; 
s[i # p] :=cos*aip+sin*aiq; 
s[i,q] :=sin*aip-cos*aiq 

end; 

akp:=fac; a[k, q] :=a[q,k] :=0 # 0; 

cs: =cos*sin; fac:=(cs+cs) *apq; 

cos:=cos*cos; sin:=sin*sin; 

a[p,p] : =co s * app+s i n* aqq+f ac ; 

a [ q , q ] : =s i n* app+cos *aqq- f ac ; 

atp* ql :=a[q,p] :=(app-aqq) *cs-apq*< cos-sin) 

end 



end ; 

a[k,p]:=:aCp,k]:=:akp 



end 
end givens; 
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TRIDIAGONALIZE SYMMETRIC MATRIX GIVSNS ttETEOD 



procedure vecgivens (a,s ,n) ; value n; integer n; array a,s ; 

comment HUCC LIBRARY PROCEDURE LE32: 
AUTHOR : J. BOOTHROYDs 

As for LE31 except that only the upper half of the 
symmetric matrix is stored, by columns , in the vector 
a[l;n*(n+l) div2 ]« The continued product V^V^Vy... 
of orthogonal transformation matrices is available in 
s[l:n,l:n]. If this procedure is used as an alternative 
to LL09 and prior to LL10 and LL11 it will be necessary 
to transfer the diagonal and codi agonal elements of array 
a to vectors c[l:n] and b[l:n] respectively. This may be 
done by 

jlessl:=k:=l; c[l] :-a[l] ; 
for j : = 2 s tep 1 until n do 
begin k:=k+j ; c[j ] :=a[k] ; 

b[ jlessl] :=a[k-l] ; jlessl:=j 

end; 

b [n] :=0.0; 

The eigenvectors of the tri diagonal system should 

be premultiplied by s [l:n,l:n] to yield the eigenvectors 

of a; 
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procedure vecgivens(a,s,n) ; value n; integer n; array a,s; 
begin integer i,p,q,nless1 ,k,cp,cq,ci,ip,iq; integer array c[1:n]; 

real akp,app,akq,apq,aqq,aip,aiq,fac,cos,sin,cs; 

for i ; =1 step 1 until n do 

begin s[i,i] :=1 .0 

for j:=i+1 step 1 until n do s[i, j] :=s[ j,i] :=0 # 

end ; - 

p:=0; for i:=1 step 1 until n do begin c[i]:=p; p:=i+p end; 
nless1:=n-1; 

for p:=2 step 1 until nlessl do 
begin k:=p-1; cp:=c[p]; akp:=a[k+cp] ; 
for q:=p+1 step 1 until n do 
'W^. begin cq:=c[q]; akq:=a[k+cq] ; 

apq : =a[p*cq] ; app : =a[p+ cp ] ; aqq : =a[ q+cq] ; 
f ac : =s qr t ( akp* akp+akq* akq) ; 
if facfcO^O then 

begin cos : =akp/f ac ; sin:=akq/fac; 
for i:=1 step 1 until n do 
begin if i>p then 

begin ci:=c[i]; ip:=p+ci; 

iq:= if i>q then q+ci else i+cq; 
aip:=a[ip] ; aiq:=a[iq]; 
a[ip] : =cos*aip+sin*aiq; 
a[iq] :=sin*aip-cos*aiq 

end; 

aip:=s[i,p]; aiq:=s[i,q] ; 
s[i,p] :=cos*aip+sin*aiq; 
s[i,q] :=sin*aip-cos*aiq 

end i; 

akp:=fac; a[k+cq] :=0„0; 
-~ cs:=cos*sin; fac:=(cs+cs)*apq; 

cos : =cos*cos ; sin:=sin*sin ; 
a[p+cp] : =cos*app+sin*aqq+f ac ; 
a[q+cq] :=sin*app+cos*aqq-fac; 
a[p+cq] : =( app-aqq) *cs -apq*( cos -sin) 

end 

end q; 

a[k+cp] :=akp 

end p 
end vecgivens ; 
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EIGENVALUES AND EIGENVECTORS OF SYMMETRIC MATRIX 

procedure jacobi (a,s ,n,rho) ; value n,rho; integer n; real rho; 

array a,s ; 
comment HUCC LIBRARY PROCEDURE LLOls 

AUTHOR ACM 85 : 

The procedure finds all eigenvalues and eigenvectors of 
a real symmetric matrix by Jacobi' s method as described 
in "Mathematical Methods for Digital Computers" edited 
by Ralston and Wilf « The eigenvectors of a[l:n,l:n] 
are built up in s[I:n,l:n] the kth eigenvector 
occupying column k. The corresponding eigenvalue 
occupies element a[k,k] of the original matrix, 
rho is the precision tolerance for the process which 
is terminated when, for every off diagonal element 
a[i , j ] , abs(a[i, j])< rho x norml/n where norml is the 
square root of the sum of the squares of the off 
diagonal elements of a; 



procedure jacobi(a,s,n,rho) ; 
value n,rho; 

integer n; real rho; array a,s; 

begin real norml , norm2 , thr , mu , omega , sint , cost , int 1 , v 1 f v2 , v3 ; 
integer i f j,p,q,ind; switch ss:=main,main1 ; 

for i:=l step t until n do 
for j:=l step 1 until i do 

if i=j then s[i f j}:=1 # else s[i, j] :=s[ j,i] :=0; 
intl :=0; 

for i:=2 step 1 until n do 

for j:=1 step 1 until i-1 do int1:=int1+2*a[i, j]t2; 
norml :=sqrt(int1) ; norm2 : =( rho/n) *norm1 ; 
thr : =norml ; ind : =0 ; 

main: thr:=thr/n; 
mainl: for q:=2 step 1 until n do 
for p:=l step | until q-1 do 
if abs(a[p,q]) > thr then 

begin indT=1; vl:=ra[p,p]; v2:=a[p,q] ; v3:=a[q,q] ; mu:=.5*(v1-v3> ; 

omega:=if mu=0 # then -1 # else -sign(mu)*v2/sqrt(v2*v2+mu*mu> ; 
sint :=omega/sqrt(2*(l+sqrt(l-omega*omega>)) ; 
cost : =sqrt< 1-sint*sint ) ; 
for i:=1 step 1 until n do 

begin int 1 : =a[i f p]*cost-a[i , q]*sint ; 

a[i , q] : =a[i , p]*sint+ a[i , q]*cost ; a[i,p] : =int 1 ; 

int1:=s[i,p]*cost-s[i,q]*sint; 

s[i,q]:=s[i,p}*sint+s[i,q]*cost; sfi,p] ;=int1 

end ; 



a 
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for step 1 until n do 

begin a[p,i]:=a[i,p]; a[q,i] := a [i, q] 
end ; 

a[p,p] :=v1*cost*cost+v3*sint*sint-2*v2*sint*cost ; 
a[q,q] : =v1*sint*sint+v3*cost*cost+2*v2*sint*cost ; 
a[p,q]:=:a[q f p] :=(v1-v3)*sint*cost+v2*(cost*cost-sint*sint) 

end ; 

if ind=1 then 

begin ind:=Q; go to maim 
end 

else if thr > norm2 then go to main 
end jacobi ; 
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SOLUTION OF THE EIGENPRQBLEM (A~AB)x=Q 

procedure eigensclve (A,B ,n ,eigen,eps ,nonposdef) ; value n,eps ; 
integer n; real eps ; array A,B,eigen; label nonposdef ; 

comment HUCC LIBRARY PROCEDURE LL 02 
Author : J, Boothroyd, 

Solves the equation (A-XB)x=0 for symmetric A,B [l:n,l:n] 

provided one of these is positive definite. The equation 

is transformed into (C-AI)y=0 where C«L~ 1 £.(L T )~ 1 

T 

is symmetric and y=L x. If B is non positive definite 
the Choleski decomposition fails in the attempted evaluation 
of a square root of a negative value . In this event the 
original equation is rearranged as (B-AA) x=0 for which the 
eigenvalues are the reciprocals of (A-AB)x-O and the 
eigenvectors are unchanged. Failure of this second attempt 
causes exit to the label nonposdef. 

At a successful exit the eigenvalues are in eigen[l:n] and 
the vectors occupy A, This procedure uses MS 01, the revised 
issues of LE 16, LE 17, together with LE 18 and LL 01. The 
parameter eps is a precision tolerance used by LL 01 for 
details of which see the appropriate commentary; 



( 
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procedure eigensolve(A,B,ii f eigen,eps ,nonposdef) ; value h,eps; int eg gr n; real eps; array A,B,eigen; label nonposdef- 
begin integer i, j; boolean recip; arra y C[l:n,l:n]; switch s:=newtry,sv/rp,ok; " ~~ * ' 

comment temporary storage of the main diagonal of B; 
for i:= 1 step 1 until n do eigen[i] :=B[i f i] ; 
recip :=false; commen t assumes B is positive definite; 
newtry: choleski(B,n, swap) ; ~ 

comment clear upper triangle of B; 

for ±: ~ 1 step i until n do for j:= i+| step 1 until n do B[i, j]s= 0.0; 
goto ok; 

swap: jLf recip then goto nonposdef; recip: =t rue ; 

comment B was nonpositive definite. Swap A,B and try again; 
for i:= 1 step 1 until n do 
begin B[i ,i] :=A[i ,i] ; A[i,i]:= eigen[i]; 
for j:= i+1 step 1 until n do 

begin B[j,i]:=A[i, jl; A[i , j] :r=A[ j, i] :=B[i , j] end 

end ; 

goto newtry; 
ok: linv(B,n) ; comment forms L-1 in B; 

mxmult(B,A,C,n,n,n) ; comment C : —L-1 A: 

comment transpose L-1 and clear lower triangle of B; 

for i:= 1 step 1 until n do 

for j:= i+i step 1 until n do begin B[i, j] :=B[ j,i] ; B[j,i]:=0.0 end; 

mxmult(C f B,A,n r n,n) ; ' comment A:= C(Lt)~1 = L-1A<Lt)-1 ; 

jacobi(A,C,n,eps); comment ACM85. Eigenvalues on A diagonal, y vectors in c- 
for i:_ 1 step l until n do eigen[i] := if recip then 1.0/A[i,i] else A[i,i]; 
mxmult(B,C,A,n,n,n) ; comment x vectors = (Lt)-ly now in A- 
end eigensolve; 1 



LL03 to LL07 

i ii i ii 

SOLVE SYMMETRIC EI GENS YS TEH 

These procedures are taken from the LINEAR ALGEBRA Handbook Series of 
Numerische Mathematique 4 pp 354 - 376 (1962), Their author, J. II. Wilkinson 
of the National Physical Laboratory, England, is the recognised authority 
on this subject and these procedures are included In the HUCC Library as 
examples of classical algorithms in this field. 

The following program was used to test four of the five procedures and 
indicates how to use the various parameters. 




TEST HOUSEHOLDER; 

begin comment procedure declarations; 
begin integer n; read n; 

begin integer i, j,mr,ml,t,m1 ; real e, gamma, norm; 
array z, a[1 : n, 1 : n] , c,b, w[1 : n] ; 
for i:=1 step 1 until n do 
for step 1 until i do read a[i, j]; 

householder ( a, n,c,b) ; 

tridibisection2(c,b,n,-20,n,0,30, 10 -16,w > norm,m1) ; 
tridiinverseiteration(c,b,n,w,norm,m1 , 10 -8,z) ; 
backtrans format ion( a , b , z , n , ml ) ; 
print ££1?EIGENVALUES£12?? ; 
sameline; scaled(9) ; 

for i:=1 step 1 until ml do print w[i] ; 
print ££15?EIGENVECTORS£l?? ; 
for i:=1 step 1 until n do 
begin print ££12??; 

for j:=1 step 1 until n do print z[j,i],££s2?? 

end i 

end 

end 

end; 



LL03 

HOUSEHOLDER TRIDIAGONALIZATION « FULL MATRIX 
procedure householder tridiagonalization(a,n,c,b) ; 
value n; integer n; array a,b,c; 
comment HUCC LIBRARY PROCEDURE LL03: 
AUTHOR : J.H. WILKINSON: 

The symmetric matrix a[l:n,l:n] is reduced to a symmetric 

tridiagonal matrix using only the lower triangle of a. The 

diagonal elements of the tridiagonal matrix are stored in c[l:n], 

the subdiagonal elements in b[lsn] with bfn^O. 

Sufficient details of the transformation are retained in 

a and b to enable the eigenvectors of a to be formed from the 

tridiagonal eigenvectors using procedure back trans formation; 



LL04 

EIGENVALUES OF SYMMETRIC TRIDIAGONAL MATRIX 1 
procedure tridibisection l(c,b ,n,gu,go, t,gamma,w,norm,ml) ; 
value n,gu,go,t, gamma; integer n,t,ml; 
real gu,go,gamma,norm; array c,b,w; 
comment HUCC LIBRARY PROCEDURE LL04: 
AUTHOR : J.H. WILKINSON t 

c[lsn] and bfl:n] are the diagonal and sub-diagonal of 
a symmetric tridiagonal matrix. The procedure determines 
the number ml of eigenvalues lying between gu and go and 
computes these eigenvalues in w[l:n] in decreasing order 
of magnitude by the method of bisection, t is the number 
of bisection steps (30 to 35is appropriate on a 503), 
norm is the infinity norm of the matrix and gamma is the 
square of the relative machine precision (10~ 16 for a 503); 



LL05 

EIGENVALUES OF SYMMETRIC TRIDIAGONAL MATRIX 2 
procedur e tridibisection aCc^jnjejmr^jt^ainmajWjnormjml) ; 
value n,e f mr f mA,t f gamma; real e,gamma s norm; 
Integer n,mr,mA t t, ml; array c 9 b»w; 
comment HUCC LIBRARY PROCEDURE LL05: 
AUTHOR : J.H. WILKINSON: 

c[lsn] ani r in] are the diagonal and subdiagonal of a 
symmetric tridiagonal matrix. The mr .eigenvalues to the left of e 
and the taS eigenvalues to the right of e are computed and stored 
in w[lsn], t is the number of bisection steps(30to35 for a 503), 
norm is the infinity norm of the matrix and gamma is the square 
of the relative machine precision (10~ 16 for 503) ; 



LL06 



EIGENVECT OR C? TRIDIAGONAL MATRIX 

procedu re tridiinverscl teration(c,b ,n,w,norm,wl,macheps ,z) ; 
value n 9 ml,norm,machcps ; integer n 5 ml; real norm,macheps ; 
array CjbjWjjZ; 

comment HUCC LIBRARY PROCEDURE LL06 : 
AUTHOR : J.H. WILKINSONS 

c[lsn] and hflsn] are the diagonal and subdiagonal of a symmetric 

tri diagonal matrix. w[l?n] contains the ml eigenvalues from w[l] 

to w[mll. norm is the infinity norm of the matrix and macheps 

is the relative machine precision (10~ 8 on a 503). The ml 

eigenvectors are computed and stored in z[l:n,l:n], element z[i,j] 
th 

being the j component of the eigenvector corresponding to the 



eigenvalue w[i]; 

LLP 7 

EIGENVECTO RS OF SYMMETRIC MATRIX 

procedure backtransfcrmritionCa^b ,z s n,ml) ; value n,ml; integer n 5 ml; array a,b, 
comment HUCC LIBRARY PROCEDURE LL07: 
AUTHOR : JJI. WILKINSON: 

b[l:n] is the subdiagonal of a symmetric tridiagonal matrix whose 
eigenvectors occupy z[l:n,lrn] ,by rows. The eigenvectors of the 
original matrix a[l:n,l:n] are derived and overwritten on array z; 



LL03 



procedure householder tridiagonalisation(a, n)result :<c,b) ; 
value n; integer n; array a,b,c; 

begin integer j,i,k; real ai,sigma,h,bj,bigk,bi ; array q[1:n-1]; 
for i : =n step -1 until 3 do 
begin sigma: =0; 

for k:=1 step 1 until i-1 do 

sigma : =sigma+a[ i , k] *a[ i , k] ; 

ai:=a[i,i-1]; 

if ai > then bi:=-sqrt( sigma) else bi :=sqrt( sigma) ; 

b[i-1]:=bi; 

it bi*0 then 

begin h:=sigma-ai*bi ; 

a[i # i-1]:=ai-bi; 

for j:=i-1 step -1 until 1 do 
begin bj:=0; ~ 

for k : =i -1 step -1 until j do 

bj:=bj+a[k, j]*a[i,k]; 

for k:=j-1 step -1 until 1 do 

bj:=bj+atj,k]* a [i,k]; 
q[j]:=bj/h 

end j; 
bigk:=0; 

for step -1 until 1 do 

bigk:=bigk+a[i,j]* q [j] ; ~ 

bigk:=bigk/(2*h) ; 

for j:=i-1 step -1 until 1 do 

q[ J 3 : =q[ J ] -bigk*a[ i , j ] ; 

fof j:=i-1 step -1 until 1 do 

begin for k:=j step -1 until 1 do 

a[j,k]:^a[j,k]-a[i,j]*q[k]-a[i,k]*q[j]; 

end j 
end 

end i; 

for i:=n step -1 until 1 do 
c[i]:=a[i,i]; "~ 
b[1]:=a[2,1]; 
b[n]:=0 



LL04 



procedure tridibisectionl (c,b,n, gu, go, t , gamma) result : (w, norm, ml ) ; 

value n, gu, go, t, gamma; 

integer n, t ,m1 ; 

real gu, go, gamma, norm; 

array c,b,w; 

begin integer i, j,k,al ,a2 f d; 

real l,g,h, lambda, pi ,q1 ,y ; switch ss : =noeigenvalue ; 

array p[1 :n] ; 

procedure sturms sequence; 
begin p1:=0; q1:=1; al :=G; 

for i : =1 step 1 until n do 

b egin y : =<c[i] -lambda) *q1 -p[i]*p1 ; 

" pi :=qi ; qi :=y; 

if pi > = ql > then a1 : =a1 +1 

end i; 

if ql =0 and pi >0 then al : =a1 -1 

end; 

if gu>go then 

begin g:=gu; gu:=go; go:=g end ; 
norm:=abs(c[1 ]) +abs(b[1 ]) ; 
for i:=2 step 1 until n do 

begin 1 : =abs(b[i-1 ])+abs(c[i])+abs(b[i]) ; 
if l>norm then norm: =1 

end ; 

for i:=1 step 1 until n-1 do 

begin if b[i]=0 then p[i+1 ] : =gamma*norm*norm 
"else p[i+1]:=b[i]*b[i] 

end ; 
p[1]:=0j 

if gu>norm or go<-norm then 

begin ml :=0; goto noeigenvalue end ; 

lambda: =gu; 

sturms sequence; 
a2 : =a1 ; 

if q1=0 then a2:=a1+1 ; 
lambda; =go; 
sturms sequence; 
ml :=a2-a1 ; 
d:=a1 ; 

if go>norm then go:=norm; 
if gu<-norm then gu:=-norm; 
for k:=1 step 1 until ml do 
be gin d: =d+1 ; 

g:=go; h:=gu; 
for j:=1 step 1 until t do 

begin lambda : =( g+h) /2 ; 
sturms sequence; 
if al > d then h:=lambda else 
g: =1 ambda 

end j; 
w[k]:=<g+h)/2 

end k; 



9 



noeigenvalue: end; 



LL05 



procedure tridibisection 2 (c,b,n,e,mr,ml, t , gamma) result :(w, norm, ml) ; 
value n,e,mr , ml, t, gamma; 
integer n,mr,ml, t,m1 ; 
real e, gamma, norm; 
array c,b,w; 

begin integer dl ,d2,i, j,k,a1 ,d; 

real l,g,h, lambda, p1 ,q1 ,y; 
array p[1 ;n] ; 

procedure sturms sequence; 

begin p1:=0; q1;=1; a1:=0; 

for i:=l step 1 until n do 
begin y:=(c£i]-lambda)*q1-p[i]*p1 ; 

Pl:=q1; q1:=y; 

^w*/ if p1 > s ql > 

then aT:=a1+1 ; 

end i; 

if q1=0 and p1>0 then al :=a1 -1 
end ; 

norm:=abs<c[1])+abs(b[l]) ; 
for i:=2 step 1 until n do 

begin l;=abs(b[i-lIT-fabs(c[i])+abs(b[i]) ; 

if l>norm then norm:=l 

end ; 

for i:=\ step 1 until n-1 do 

begin if b[i]=0 then p[i+1 ]::sgamma*norm*norm 

else p[i+1]:=b[i]*b[i] 

end; 

pD3:=0; 
lambda: =e; 
sturms sequence; 

' W 

dl : =a1 -mr ; d2 : =a1 +ml ; 
if dKO the n dl :=0; 
if d2 > n+1 then d2:=n; 
d:=d1 ; m1:=0; 

for k:=d! step 1 until d2-1 do 
begin d:=d+1; 

g : =norm ; h : =-norm ; 

for j:=1 step 1 until t do 

begin 1 ambda:=(g+h) /2 ; 

sturms sequence; 

if al > d then h:=lambda else g:=lambda 

end j; 

ml : =m1 +1 ; 

w[m1]:=(g+h)/2 

end k 



end; 



LL06 



procedure tridiinverse iteration(c,b,n,w f normal ,macheps)results:(z) ; 
value n,m1 ,horm,macheps; integer n,m1 ; real norm,macheps; array c,b,w,z 
begin integer i,j; real bi,bi1 ,z1 .lambda, u,s,v,h,eps, eta; 
array m,p, q,r,int[1 :n] ,x[l :n+2] ; 
1 ambda : =norm ; eps : =macheps*norm ; 
for j:=1 step 1 until ml do 
begin lambdas =lambda~sps ; 

if w[j]<l ambda then lambda:=w[ j] ; 
u: =c[1 ]-lambda; v:=b[1 ] ; 
if v=0 then v:=eps; 
for i:=1 step 1 until n-1 do 
begin bi:=b[i]; 

if bi=0 then bi:=eps; 
bil :=b[i+1]; 
if bi1=0 then bil :=eps; 
if >Jbs( bi> > abs(u) then 
^egin m[i+ll:= u /bi; 

if m[i+1]=0 and abs(bi) < eps then m[i+1]:=1 ; 

p[i] :=bi ; q[i] :=c[i+1 ]-l£mbda; 

r[i]:=bi1 ; 

U2=v-m[i+1]*q[i] ; 

v:=-m[i+1]*r[i]; 

int[i+1]:=1 

end 

els® begin m[i+1 ] :=bi/u; 
p[i]:=u; q[i]:=v; 
r[i]:=0; 

u:=o[i+1 3-lambda-m[i+1 ]*v; 
v:=bi1 ; int[i+1 ] :=-l 

end 

end i; 

p[n] : =u ; q[n] : =r[n] : =0 ; 

x[h+1 ] : rrxCn+2] t =0 ; h : =0 ; et a : =1 /n ; 

for i : =n £tep -1 until 1 do 

^ £21*2 ^=®*a-qCi]*x[i+1]-r[i]*x[i+2]; 

£f p[i]=0 then x[i]:=u/eps 

else x[i]:=u/p[i]; 
h:=h+abs{x[i]) 

end i ; 
h:=1/h; 

for i:=1 £tep i until n do 

x[i]:=x[i]*h7 ~~ 

for i : =2 step 1 until n do 

begin if int[i]>0 then 

£3g£n u:=3:[i-1]; 

x[i-1]:=x[i]; 

x[i]:=u-m[i]*x[i-l] 

end 

else x[i]:=x[i]~m[i]*x[i-1] 

end i ; 

h:=0; 

for i : =n jstep -1 until 1 do 

begin ^-=xCi]-q[i]*x[i+1]-r[i]*x[i+2]; 
if p[i]=0 then x[i]s=u/eps 

else x[i]:=u/p[i]; 
h:=h+x[i]*x[i] 

end i; 
h:=1/sqrt(h> ; 

for i:=1 ste p 1 until n do 
z[j,i]:=x[i3*h "~ 

end j 



procedure backtransf ormat ion( a ? b , z , n, ml ) ; 
value h,m1 ; integer n,m1 ; array a,b,z; 

begin integer i, j,k; real s; for j:=1 step i until ml do 
for k:=3 step 1 unti l n do 
if ■b[k-1]+0 then 
begin s:=0; • 

for i:=1 step 1 until k-1 do 

s:=s+a[k,i]*z[j,i] ; 

s:=^/<b[k-1]* a [k,k-l]); , 

for i;=1 step 1 until k-1 do 

z[ j,i] :=z[ j ; i]+s*a[k,i] 

end ' 

end ; 



LL08 



EIGENVALUES AND VECTORS OF SYMME TRI C ?4ATRI X 

procedure vecjacobi(a,s ,n,rho) ; value n,rho; real rho; 
integer n; array a,s ; 

comment I1UCC LIBRARY PROCEDURE LL08; 
AUTHOR : J. BOOTHROYD: 

A revised version of LL01, jacobi , such that only the 
upper half of the symmetric matrix is stored in the 
vector a[l:n*(n+l)div2] . Storage is by columns as 
described in LE24. On exit from the procedure the 
eigenvalues occupy elements a[l] through a[n] and the 
eigenvectors occupy the array s[l:n,l:n] . On a 20 x 20 
matrix LL01 took 120 seconds, LL08 only 47 seconds. On 
the same matrix J. II. Wilkinson's procedures (LL03,04,06,07) 
took 27 seconds and their improved counterparts (LL09, 10,11, 
12) took 23 seconds; 



LL08 



procedure vecjacobi(a,s,n,rho) ; value n,rho; real rho; integer n; array a,s; 
begin integer array c[l:n]; integer i , j , ci , c j,p, q,cp,cq, jlessl , qlessl , ip, iq; 
switch ss :=main,main1 ; 

real fac,ai j, thr,norm1 ,norm2, apq, app, aqq,m,mu, lambda, cost, sint,aip,aiq, sip, siq,sincos; 

boolean ind; 

p : =0 ; fac:=0,0; 

for i : =1 step 1 un til n do 

begin s[i ,i j :=1 c 0; c[i]:=p; p—p+i; 

for j:=i+1 step 1 until n do s[i, j] :=s[j,i]:=0..0; 

end ; 

for j:=2 £tep 1 unti l n do 
begin cj:=c[j] ; jlessl :=j-1 ; 

for i : =1 stop 1 unt il jlessl do 

begin ai j : =a[i+c j] ; f ac : =f ac+2, 0*ai j*ai j end 

thr : =norm1 ; =sqrt( f ac) ; nor:n2 : =rho*norm1 /n ; 
main: thr :=thr/n; 
mainl : ind:=f alse; 

for q:=2 step 1 until n do 
begin cq: =c[q] ; qleasl :=q-1 ; 

for p:=1 step 1 until qlessl do 
begin apq:=a[p+cq] ; 

if abs(apq) > thr the n 
begin cp:=c[p] ; ind: = true ; 

app:=:a[p:-cp]; aqq:=a[q+cq] ; m:=app-aqq; 
lanbda:=sign(m) *apq; mu:=0.5*abs(m) ; 
fao:=Q t) 5/sqrt(lambda*lambda+mu*mu) ; 
cost :=sqrt(0 . 5+nu*fac) ; sint : =lambda*f ac/cost ; 
for i:=1 step 1 until n do 
begin ci : =c[i] ; 

if i < P then begin ip:=i+cp; iq:=i+cq end 
else begin ip:=p+ci; 

iq:=if i>q then q+ci else i-fcq 

end ; 

aip:=a[ip] ; aiq:=a[iq]; 
sip:=s[i,p] ; siq:=s[i,q] ; 
s[i,p] :=cost*sip+sint*siq; 
s[i , q] : =sint*sip-cost*siq; 
a[ip] : =cost*aip+sint*aiq; 
a[iq] : =sint*aip-cost*aiq 

end i; 

sincos :=sint*cost ; f ac : =( apq+apq) *sincos ; 
sint :=sint*sint ; cost : =cost*cost ; 
aCp+cp] : =cost*app+sint*aqq+f ac; 
a[q+cq] : =sint*app+cost*aqq-f ac ; 
aCp+cq] :~0 o 
end 

end 

end; 

if ind then goto mainl else if thr>norm2 then goto main; 
for i:=2 step 1 until n do a[i] :=a[i+c[i] ] 
end vecjacobi; 
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£££codurj householder (a, n,c,b); value n; integer n; real array a, b,c; 

38368 HH r ^ j ^^ 32>iPlU3l>Cj>Ck>kCJ>iplUSl ' nle ^ T ^aij^i.h^ho^j; 
nless2:=n-2; nlessl :=n-1 ; 
tor i:=1 step 1 until nless2 do 

begin iplus1:=i+l; s:=0,0; k:'=ciplus1 : =col[ipl U 8l ] ; 
for j:=iplus1 step 1 until n do 
begin aij:=a[i+k]7 k:=fc+j; s:=s+ai j*ai j end; 
ai j;=a[i+ciplus1 ] ; bi:= £qr t(s) ; — m 
if aij > 0,0 then bi:=-bi; bFi]:=bi; 
if bi £ 0.0 then 

begin h:=Q-aij*bi; a[i+ciplus1 ]:=ai j-bi; 
£££ j:=iplus1 step 1 until n do 
bejin s:=0 8 0; cj:=coT[jT7 ~ 

for k:=iplus1 step 1 until n do 
begin ck:=col[k]; " 

s:=S4-a[i+ck]*a[if k<j then k+cj else j+ck] 
end; "~ — — — 

z[ jj :=s/h 

end; 

s:=0 # 0; k:=ciplus1 ; 

l^.i/^ 1 ^ 1 ^ 1 n & b££in s:=s + a[i + k]*z[j]; k:=k + j end; 
rho:=s/(h+h) ; k:=ciplus1 ; 9 

|££ J|=ipluBl step 1 until n do begin «[j]:«[j]^ho*ati+k]j k:=k + j end; 

±21 J2=iplus1 step i until n do 

£SSi2 cj:=colCj]; ai j:==a[i +C j] ; zj:= z [j] ; ck:=ciplus1; 

£or k:=iplus1 stop 1 until j do 

begin kcj:=k+ c j; m 

a[kcj]:=a[kcj]-a[i+ck]*zj-aij*2[k]; ck:=ck+k 

end 

end 

c[i] :=a[i+col[i]] 

_end ; 

cj:=col[n]; c[nle3sl]:= a [nless1+col[nless1]]- 
b[nles8l]:=aCnl«ssl+cj]; b[nl:=:0.0; c[n] : =a[n+c j] 
end householder; JJ 



LL09 to LL12 

SYMMETRIC EIGENSYSTEM HALF MATRIX 

These procedures are adapted from LL03 to LL07 to compute the eigenvalues 

th 

and eigenvectors of an n order symmetric matrix of which only the upper 
half is stored by columns in a vector a[lsn*(n+l) div2 1« In designing 
this revision the opportunity has been taken to improve, where possible, 
the efficiency of the original procedures and rationalise, to some extent, 
the parameter organisation. 

The following program was used to test LL09 to LL12 and illustrates the 
use of the several parameters 



TEST HOUSEHOLDER ; 

begin integer n,t; real eps; switch ss:=L; 
L: read n,eps,t; 

begin integer i, j,k; real norm; 

array z[1 :n,1 :n],a[1 :n*(n+1) div 2],c f b,w[1 :n]; integer array col[1:n]; 

comment procedure declarations; 

k:=0; 

for j:=1 step 1 until n do 
begin col[j]:=k; 

for i:=1 step 1 until j do read a[i+k]; 

k:=k+j 

end ; 

householder(a,n,c,b) ; 
eigbisec(c,b,n,eps,t,w,norm) ; 
trivectors(c,b,n,w,norm,eps,z) ; 
eigvectors(a,b,z,n) ; 
freepbint(9) ; sameline; 
for j:=1 step 1 until n do 
begin print ££!??, w[j] ,££!??; 

for i : =1 step 1 until n do 

begin .if i-i div 10*10=1 then print ££!??; print z[i, j] end 



end 



end ; 
goto 



LL09 

HOUSEHOLDER TRIDIAGONALISATION - HALF MATRIX 

procedure householder(a.n T c.b'> j value n; integer n; array a,c,b; 

comment HUCC LIBRARY PROCEDURE LL09; 

AUTHOR : J, BOOTHROYD: 
th 

Transforms an n order symmetric matrix whose upper half is stored 
by columns in a[l:n*(n+l) div3j to tridiagonal form he 
diagonal and super-diagonal Stored in c[lm] and b[l:n] 
respectively (with b[n]=0). Sufficient information is retained in 
a, and b for subsequent use by procedure eigvectors (LL12) ; 

LL10 

EIGENVALUES OF TRIDIAGONAL MATRIX 

procedure eigbisec(c,b,n,eps,t,w,norm) ; value n»eT>s.t: 
integer n,t; real eps,norm; array c,b,w; 
comment HUCC LIBRARY PROCEDURE LL10: 
AUTHOR : J. BOOTHROYD: 

Evaluates, in decreasing order of magnitude, using the method of 
bisection, the n eigenvectors of a tridiagonal matrix whose 
diagonal and co-diagonal occupy c[l:n] and b[l:n] respectively. 
The results occupy w[l:n]. eps is the relative machine precision 
(10~ 8 on a 503), t is the number of bisection steps (30 to 35 is 
appropriate) and norm is the computed infinity norm of the matrix; 

LL11 

EIGENVECTORS OF TRIDIAGONAL MATRIX 

procedure trivector(c,b ,n,x<r, norm, eps ,z) ; value n,eps,norm; 
integer n; real norm, eps; array c,b,w,z; 
comment HUCC LIBRARY PROCEDURE LLllj 
AUTHOR s J, 300THR0YD: 

Computes in the columns of z[l:n,l:n] the eigenvectors of a symmetric 
tridiagonal matrix whose eigenvalues occupy w[l:n] in decreasing 
order of magnitude f c[l;n] and b[l:n] are the diagonal and co- 
diagonal of the tridiagonal matrix* eps is the relative machine 
precision (10~ 8 on a 503), norm is the infinity norm of the matrix; 



LL12 

EIGENVECTORS OF SYMMETRIC MATRIX (UPPER HALF) 
procedure eigvectors(a,b ,z,n) ; value n; 
integ er n; array a,b,z; 
conraent HUCC LIBRARY PROCEDURE LL12; 
AUTHOR : J a BOOTHROYD: 

From the information retained in a[l:n*(n+l)div2] and b[lsn] 
following the use of LL09 and the eigenvectors of the tridiagonal 
matrix in z[l:n 9 l:n] 9 computes the eigenvectors of the original 
matrix aj 




»" 1 r ll f 
LL10 



procedure eigbisec(c, b,n, ops, t,w, norm) ; value n,eps,t; integer n,t; 
real ops, norm; real arra y b,c f w; 

begin integer i,k,a,j; real g,h,p1 ,q1 , bi,lim, lambda, y; array p[l:n]; 
g:=abs(b[1]) ; pi :=abs(c[1 ])+g; 
for i:=2 step 1 until n do 
begin h:=abs(b[i]) ; ql :=g+h+abs(c[i]) ; 
if q1>p1 thon p1:=q1; g:=h 

end ; 

lambda: =norm:=lim;=p1 ; bi:=b[1]; 
for i : =2 £tep 1 until n do 
begin if bi=0,0 then bi:=eps*lim; 
p[i]:=bi*bi; bi:=b[i] 

end ; 

for k:=1 £tep 1 until n do 
begin g:=Tambda; h:=-lim; 

for step 1 until t do 

begin lambda: =(gn3iT72 # 0; 

p1:=0 e 0; q1:=1„0; a:=0; 
for i : =1 step 1 until n do 
begin y s =( c [ i ] -1 ambda) *q1 -p[ i ] *p1 ; 
P"l:=q1; q1:=y; 

if p1>0 o = ql > f then a:=a+1 
end i ; ~" 

if q1=0 t and p1>0 then a:=a-1 ; 
if a>k then h;=l ambda else g:=lambda; 
if g=h thon j:=t 

end j; 

w[k]:=(g+h)/2 e 

end 

end eigbisec; 



LL11 



procedure trivectors(c,b,n,w,norm,eps,z) ; value n,eps,norm; 
integer n; real norm,eps; array c,b,w,z; 

begin integer i, j, iplusl ,iless1 f nless1 ; array m,p,q,r,int,x[1 :n] ; 
real bi,bi1 , lambda, u,v,xi, miplusl ,h,eta; 
lambda: =norm; eps:r=eps*norm; nlessl :=n-1 ; 
q[n]:=r[n]:=0,0; eta:=1 *Q/n; 
for step 1 until n do 

begin lambda: =lambda~eps; 

if w[j]<lambda then lambda:=w[ j] ; 
u:=c[1 ]-lambda; bi:=*r:=b[1 ]; 
if bi=0 o then bi:=v:=eps; 
for i:=1 step 1 until nlessl do 
begin iplusl : =i+1 ; 
^ bi1:=b[iplusl]; if bil=0.0 then bi1:=eps; 

if abs(bi)>abs(u) then 

begin miplusl : =m[ iplusl ] :=if u=0.0 and abs<bi)<eps then 1 # else u/bi; 
p[i]:=bi; bi:=r[i]:=bi1 ; 
bil : =q[i] :=c[iplus1 ]-lambda; 
u : =v-miplus1 *bi1 ; v : =-miplus1 *bi ; 
int[iplus1 ] :=1 «0 

end 
else 

begin miplusl : =m[ iplusl ] :=bi/u; 

p[i]:=u; q[i]:=v; r[i]:=0.0; 
u :=c[ iplusl ]-lambda-miplus1 *v; 
bi:=v:=bi1 ; int [ iplusl ] :=-1 .0 

end 

end i; 

p[n]:=if u4:0 e then u else eps; 
v : =u : =h : =0 o ; 

for i:=n step -1 until 1 do 
begin xi:=x[i] s=(eta-q[i]*u-r[i]*v)/p[i] ; 
^-^ v:=u; u:=xi; h:=h+abs(xi) 

end ; 

* — - i:=1 ££f£ 1 until n j*2 x[i]:=x[i]/h; 
u:=x[1]; ileseltsl; 
for i : =2 step 1 until n do 

begin if int[i]>0.0 then begin v:=x[iless1 ] :=x[i] ; 

u:=x[i]:=u-m[i]*v 

end 

else u:=x[i]:=x[i]-m[i]*u; 

iless1:=i 

end i; 

v:=u:=h:=0.0; 

for i:=n step -1 until 1 do 
begin xi:=x[i]:=(x[i]-q[i]*u-r[i]*v)/p[i]; 
v:=u; u:=xi; h:=h+xi*xi; 

end ; 

h:=T,0/sqrt(h) ; 

for. i:=l step 1 until n do z[i, j]:=x[i]*h 

end j 
end trivectors; 



c 
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procedure eigvectors(a,b,z,n) ; value n; integer n; array a # b,z; 
begin integer i,j,k,iplus1 ,ck,ciplusl ; real s,bi; 
for j:=1 step 1 until n do 
begin for i : =n-2 step -1 until 1 do 
begin bi:=b[i]; 

if bi+0.0 then 

begin s:=0.0; iplus1:=i+1; ck:=ciplusl :=col[iplusl ] ; 

k:=iplusl step 1 until n do begin s:=s+z[k, j]*a[i+ck] ; ck:=ck+k end ; 
s : =s/(bi*a[i+ciplus1 ]) ; ck: =ciplus1 ; 

for k:=iplusl step 1 until n do begin z[k, j]:=z[k, j]+a[i+ck]* s; ck:=ck+k end 

end 

end 

end 

end eigvectors; 



LZ 01 



GENERATE UNSYMMETRIC TEST MATRIX 

procedures testmx(a,n) ; value n; Integer n; array a; 

comment HUCC LIBRARY PROCEDURE LZ 01: 
AUTHOR J, B00THR0YD: 



The procedure generates matrices of the type described by T, J.Dekker, 
Report No. MR 63, Mathematical Centre, Amsterdam. 
The matrices have the following properties 

(a) elements a[i , j ] are integers 

(b) elements of the inverse, (-1) i+ ^*a[i , j ] 9 are also integers 

(c) the degree of ill condition increases rapidly with n 

(d) the determinant of all matrices is 1. 

Computation of a matrix order 15 is possible on the 503 with a 
39 bit register. As real matrices however the order of the largest 
useable matrix should be restricted to that value of n(12) for which 
all elements of the matrix have an exact floating point representation 



LZ 01 

procedure tee tmx(a,n) ; value n; integer n; arrag; a; 

begin i nteger i , j ,k, fi ,gi ,d,q,r; boolean even; integer arra^ f ,g[1 :n]; 
cemmerjt first we compute F = diag(fi) ; 
fi:=f [T] :=n; j:=n*n; 
£uj? i:=1 s_top 1 until n-1 do 
begin d:=i*i; k:=j-d; 

q:=fi div d; r:=f i-q*d; 
f [i+1 ] :=fi:=q*k+(r*k) div d 

end; 

co men t and now, using a modified prime factors algorithm to obtain G = diag(gi) we compute 
FG-1 , whose elements replace those of F; 
icr i:=1 step 1 until n do 

begin d: = gi:z:1 ; qs=f i :=f [i] ; j:=2; 
"icwj: even:= jf a 1 s e ; 

next: if q > j then 

begin q:~f i div j; 

if fi4=q*j ^then begin j := j+d ; d:=2; goto newj end ; 
jif even then gi:=gi*j; evenir not even; 
f i :=q; goto next 

end; 

g[i]:=gi; f[i]:=f[i] div gi 

end; 

cooment finally , in one operation (FG-I)ffj where H is a non-existent Hilbert 

matrix whose reciprocal elements ,i+j-1 ,are computed as we go; 
for i:~1 steg 1 until n d£ 
begin fi:=f [i] ; ^ 

for j:=1 step 1 until n do 

begin gi:=g[j]; k:=i+j-1 ; 

q:=fi div k; r;=fi-q#k; 
a[i , j] : = q*gi+(r*gi) div k 

ond 

end 

GXKl tostmx; 



MC01 



LEAST SQUARES POLYNOMIAL FIT 



procedure LSQFIT(x,y ,N,a,n) ; value N,n; integer N,n; 
real array x,y,a; 

comment HUCC LIBRARY PROCEDURE MCOl: 
AUTHOR J. BOOTHROYD t 



Given the N+l sample pairs (x ,y ) ,(xi ,Yi) .... (x N ,y N ) 
in arrays x,y [0:N] the procedure determines the 

coefficients »a2» a 3 a n+l °^ t ^ e nt ^ orc * er 

polynomial approximation 

y = ai + a 2 x + a 3 x 2 + ^i+l^ 

obtained from applying the least squares principle to the 

given data. The coefficients a^a? a n+l occu Py 

the corresponding element positions of a [ 1 : n+l ] . The 
procedure uses SOLVE (LE02) to obtain the solutions of the 
normal equations . The matrix of coefficients of the 
normal equations becomes increasingly ill-conditioned 
as n increases, and it is recommended that the use of 
this procedure be restricted to values of n < 6; 



procedure LSQFIT(x,y,N,a,n) ; value N,n; integer N,n; real array x,y,a; 
be S in integer i, j,k,nplue1,nplus2; real p,q,xk; real array b[i :n+l ,1 :n+2] ; 

nplusi :=n+l ; nplus2:=n+2; 
for i:=l step 1 until nplusi do 
begin b[i,nplus2] :=0.0; 

for step 1 until i do b[i,j];=0.0 

end ; 

for k:=0 step 1 until N do 

begin p:=1 # 0; xk:=x[k]; 

to? i:=T step 1 until nplust do 
begin q:=1.0; 

for step 1 until i do 

begin b[i, j):=b[i, j]+p*q; q:=q*xk end ; 
bf i , nplus2] : =b[i , nplus2]+p*y[k] ; 
p:=p*xk 

end 

end ; 

for i:=2 step i until nplusi do 

for j:=i-l step -1 until 1 do b[j,i]:=b[i, j] ; 
SOLVEQ(b f nplusi) ; 
for i:=l step 1 until nplusi do a[i] :=b[i f nplus2] 
end LSQFIT; — 



I 



MCQ1 



HC02 

Least Squares Linear Fit 

procedure Linf it(N,x,y,aO,al,xm,ym) ; value N; integer N; 
real aO,al,xm,ym; real array x,y ; 
comment HUCC LIBRARY PROCEDURE MC02 
AUTHOR : J .11. BAXTER 

Given the N+l sample pairs (x ,y ) , (*\ 9 y\) 

in arrays x,y [0:N] the procedure determines the coefficients 
aO and al of the straight line approximation 
y=aO+alx 

obtained from applying the least squares principle to 
the given data. This is simpler than using procedure 
LSQFIT(MCOl) with n=l for this purpose. 

The procedure also evaluates the arithmetic means, 
sm and ym, of the contents of the arrays x and y 
respectively, for the reason that the straight line of 
best fit may be expressed conveniently as 

(y-ym) =al (x-xm) 
Should this not be desired, and the values of xm and ym 
be of no interest, the procedure call can be such as 

Linfit(N, independent variable, dependent variable, 

intercept, slope, slope, slope) 
so that the values of xm and ym in turn are assigned to the 
variable slope and then overwritten by the value of al; 



MC02 



procedure Linfit(N,x f y,aC ,a1 ,xm,ym); value N; integer N; real aU,al ,xm,ym; real array x f y; 
begin real n,xi,yi f sx,sy,sxx,sxy; integer i; 
n:^N+1.0; sx:=sy :=sxx::=sxy:=0.0; 
for i:=0 step 1 until N do 

begin xi:=x£i]; yi:=y[i]; 

sx:=sx+xi; sy:=sy+yi; sxx:=sxx+xi*xi; sxy :=sxy+xi*yi 

end ; 

xm:=sx/n; ym:=sy/n; 

al :=(n*sxy-sx*sy)/<n*sxx-sx*sx) ; a0:=(sy-a1*sx)/n 
end Linfit; 



) 

t 



MC03 

Weighted Linear Least Squares Fit 

procedure wt linfit(N,xk,yk,k,wk,aO,al,xm,ym) ; value N; integer N.k: 
real xk,yk,wk,aO,al,xm,ym; 
comment HUCC LIBRARY PROCEDURE MC03 
AUTHOR : J.N. BAXTER 

In cases where a curve is to be fitted to sample pairs 
(x^y^ by means of a transformation into (f(x),f f (y)) 
and a linear fit on the functions, it is necessary to 
weight the data so that the least squares approximation 
obtain is that for the original x and y. This is a Jensen 
procedure for the purpose, and k is the Jensen parameter 
which is utilised in the procedure call. For example, 
if f(x) » l /x and f f (y) » Jin y with a chosen weighting 
for each point of y*y the procedure call might be 
wtlinfit(N,1.0/x[k],£n(y[k]),k,y[k]+2) 
to find: (cept, slope, meanfx,meanfy) . 
The procedure is also capable of an unweighted fit, if 
the formal parameter wk is replaced in the procedure 
call by 1.0, but an unweighted fit is more simply carried 
out by procedure Linfit(MC02) ; 



MC03 



procedure wtlinf it (N,xk,yk,k,wk,aC ,a1 ,xm,ym) ; value N; integer N,k; real xk,yk,wk # aC,a1 7 xm,ym; 
begin real w,x,yw,xw, sw, sxw, syvv, sxxw, sxyw; , 
sxw : =syw : ^sxxw : =sxyw : =sw :=0 . ; 
for k:=0 step 1 until N do 



begin x:=xk; w : ~wk ; yw:=yk*w; xw:=x*w; 

sxw : — sxw-fxw ; sy\v:~syw+yw; sw:=sw+w; sxxw:=sxxw+x*xw; sxyw ; =sxyw+yw*x 

end ; 

xm; -sxw/sw; ym : —syw/sw ; 

al :=(sw*sxyw-sxw*syw)/(sw*sxxw-sxw*sxw) ; a0:= (syw-a1*sxw)/sw 



end wtlinf it ; 





MC04 

Standard Deviation of Fitted Straight Line 

real procedure Lindev(N,x,y ,a0 ,al) ; value N,aO,al; integer N; 
real aO,al; real array x,y; 
comment HUCC LIBRARY PROCEDURE MC04 , 
AUTHOR : J.N. BAXTER 

A procedure to evaluate the standard deviation of the N+l 
data points held in x,y[0sn] from the fitted straight line 
y~aO+alx. (The values of aO and al supplied to this procedure 
may be those calculated by Linfit(MC02) or may be other 
values, chosen at will, or derived from other work). 
To avoid the infinite value which would otherwise result 
from the use of this procedure with 2(N=1) data points 
the assignment Lindevz^S^lO 75 occurs. The procedure is 
not protected from use with N=0 which case fails on a 
sqrt error; 



real procedure Lindev(N,x,y,aO,a1) ; value N,a0,a1; integer N; real a ,a1; real array x,y 
begin integer i,n2; real sum,dev; 
sum:=0.0; 

for i:= step 1 until N do begin dev:=aQ+a1*x[i]-y[i] ; sum:=sum+dev*dev end ; 
n2:=N-1; . . 

Lindev:= (if n2±0 then sqrt(sum/n2) else 5.0 10 75) 
end Lindev; 



MD01 



LOCATE MINIMUM OF FUNCTION f (x) 



real procedure MINX(a«b .eps <,x % fx«fval) : value a,b,eps ; 
real a,b ,eps ,x,fx,fval ; 
comment HUCC LIBRARY PROCEDURE MDOl: 
AUTHOR J. BOOTHROYD ; 



A procedure to locate the minimum of fx in the interval 
a <. x ^ b by the method of trisection. fx must be 
monotonic decreasing from x=a to the position of the 
minimum and thereafter monotonic increasing until x=b. 
eps is an absolute tolerance and the procedure aims to 
locate the minimum with an error less than eps. Whether 
it succeeds will depend on how well the minimum is 
defined taking into account the specified eps and the 
precision to which fx may be evaluated for any given real 
number representation, A relative tolerance delta may 
be specified using the call MINX(a,b, (b-a)*delta,x,fx,fval) . 
On exit from the procedure fval yields the minimum value 
of the function; 



real procedure MI NX< a , b , eps , x , f x , f val ) ; value a,b,eps; real a, b, eps , x, f x, f val 
begin real sep,fx1,fx2; integer d; switch s:=L, movea; 

L: sep:=(b-a)/3.0; x:=a+sep; fxl:=fx; x:=b-sep; fx2:=fx; 
if fx2=fx1 then begin d : =0 ; b:=b-sep; goto movea end ; 
d:=i j 

if fx2>fxi then b:=b-sep else movea: a:=a+sep; 
if sep>eps then goto L; 

MINX:=x:= if d£Q then a+sep else (a+b)/2.0; fval:= fx 

end MINX; 
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CONVERT SEXAGESIMAL ANGLES TO RAD IANS 

real procedure sexrad(deg,min, sec) ; value deg,min,sec; 
integer deg,min; real sec; 

comment HUCC LIBRARY PROCEDURE MG 01: 
AUTHOR J, BOOTHROYD: 

Converts angular measure in the sexagesimal (360,60,50) system 
tc radians; 



MG 02 



CONVE RT RADIAN ANGULAR ME ASURE TO SEXAGESIMAL UNITS 

procedure radtosex(rad,deg,min,sec) ; value rad; 
integer deg,min; re al rad, sec; 

c omment, HUCC LIBRARY PROCEDURE MG 02: 
AUTHOR J, BOOTHROYD: 

Converts angular measure in radians to degrees , minutes and 
seconds in the sexagesimal system. The result is correctly 
rounded to two decimal places of seconds ; 



MG 03 



CONVERT CENTESIMAL ANGLES TO RADIANS 

real procedure centrad (deg ,min,sec) ; value deg,min,sec; 
integer deg,min; real sec ; 

comment HUCC LIBRARY PROCEDURE MG 03: 
AUTHOR J, BOOTHROYD: 

Converts angular measure from the centesimal (400,100,100) system 
to radians; 



MG 04 
■ .1. 1,,, ...... i- 

CONVERT RADIAN ANGULAR MEASURE TO CENTESIMAL U NITS 

procedure radtocent(rad , cleg, min, sec) ; value rad; 
integer deg, min; real rad,sec; 

comment HUCC LIBRARY PROCEDURE MG 04: 
AUTHOR J. BOOTHROYD: 

Converts angular measure from radians to centesimal (400,100,100) 
units. The result is correctly rounded to two decimal places of seconds ; 



MG 05 



COMPUTE DISTANCE AND GRID BEARING 

procedure bearing(x,y ? r, t) ; value x,y; real x,y,r, t; 

comment HUCC LIBRARY PROCEDURE MG 05: 
AUTHOR J, BOOTHROYD: 

Computes r, the length, and t , the clockwise radian angular bearing 
from, north of a line joining the origin (0,0) to the point (x,y) . 
For the bearing and length of a line joining (xl,yl) to (x2 ,y2) use 
the call bearing (x2~xl ,y2~yl , r , t) : 



MG 06 



COMPUTE POLAR COORDINAT ES OF PO INT X,Y 



procedure polar (x,y,r, t) ; value x,y; real x ,y,r,t; 

comment HUCC LIBRARY PROCEDURE MG 06: 

AUTHOR NATIONAL PHYSICAL LABORATORY : 



Computes the polar coordinates r,t corresponding to given x,y 
Cartesian coordinates; 



MG 07 



COMPUTE POINT OF INTERSECTION FROM ANGLES 

procedure intang (xl ,yl , a ,x2 ,y2 5 b ,x ,y) ; value xl ,yl , a,x2 ,y2 ,b ; 
real xl ,yl , a,x2 ,y2 ,b ,x,y ; 

comment HUCC LIBRARY PROCEDURE MG 07: 
AUTHOR J. BOOTHROYD: 

a 

y 

Competes the coordinates cf the point of intersection of the lines 
through (xl,yl) , (x2,y2) which make angles , with respect to the line 




(xl^yl) 

The following convention of signs nius_t be observed. 

Positive angles a are those swept out by closing the line IP onto 

12 c lockwise . 

Positive angles b are those sxv^ept out by closing the line 21 onto 
2P clockwise. 



The formulae used are : 



((x2~xl) cot(a)-(y2-vl)) 
cot(a)+cot(b) 



y= 



((y2-yl)cot(a) + (x2-xl)) 
cot(a)+cot(b) 



+yl 



< 



real proceaure soxrad(dog f min,soc) ; value deg ,min,soc ; 
ixxtog_or Jog,min; real ogc; 

bogin soc:_ ((sec/60 # 0+min) /60 % 0+abs(dog) ) *0 .01 74532925 ; 

Goxrac:=if dog<0 then -sec elso soc 
ond soxrad; 



pr ocedure raatos€x(rad,deg,min,sec) ; value rad; 
j ntoger deg,min; real pad, sec; 
begin real x; 

= x:_ ~abo(rad)/0 .01 74532925+0 .000001 39; 

d3p:= e.itior(x) ; x:= (x-deg) *6C .0 ; 

if rad<° # then deg -deg; 

min; - entier(x) ; sec:=(x~min) *60 o 0*100 .0 ; 
SGc;=entier(£sec) /1 00 .0 
oiid radtosox; 



real prcc^cure contxad(deg,min>sec) ; value dog ,min,sec ; 
i.rrt og_er ceg,mxn; real see; 

begin sec: = (<3oc/10( ,0+min)/100 B 0+abs(dog)) *0 c 01 570796327 

centrad:- i^f deg<0 then -sec else sec 
end contrad; 



procojure rau t ocoi-t ( rad , deg , mi n , s ec) ; value rad; 
integer ueg,n5n; real rad, sec ; 
bopLii real X ; 

x:= abs(rad) /0 c 01 570796327+0 o 0000005 ; 

deg:- entiex(x); x:=(x-deg) #100 o ; 

if rad<0 «0 then deg:=:«deg; 

Ldn: = entiex(x) ; sec:=(x-zain) *1 00 o 0*1 00 ; 
sec ;=oncier ( sec) /1 00 .0 
ond radtocent; 



» 
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yocec.ur e bearing(x,y ,r , t) ; value x,y ; real x,y,r,t ; 

begin real pi,piby2; pi:=3 .1 41592654; piby2:=1 .570796327; 
r :=:8qrt(x*x+y*y) ; 
if r=0.0 then t;=0.0 

else begin t:= i_f abs(x)<abs(y) then arctan(x/y) 

else piby 2-ar ctan( y/x) ; 

if y<0 .0 or x<0 .0 then t: = t+pi; 
if yX) ,0 and x<0 # then t: = t+pi 

end 

end bearing; 



m 05 



pro co du re p jlar(x ,y ,r , t) ; value x,y ; real x,y,r,t; 

beg^r real pi, piby 2; pi: =3 .1 41592654; piby 2: =1 .570796327; 
r:=sqrt(x*x+y*y) ; 
M r=0.0 then t:=0.0 

else begin t :- if abs(x)<abs(y) then piby2-arctan(x/y) 

else arctan(y/x) ; 

x:=x+y; 

if x < 0.0 or y<0 the n. t:=t-fpi; 
^f x > 0.0 and y<0 then t:=t+pi 

end 

eud polar; 



m oe 



procedaro intang(x1 ,yl ,a,x2,y2,b,x, y) ; value xl ,y1 ,a,x2,y2,b; 

real #1 ,y1 ,x2,y2,a,b,x,y ; 

^begin a:^l .0/tan(a) ; b:=a+1 o 0/tan(b) ; 

X2:=x2-x1 ; y2:=y2-y1 ; 

x:=(x2*e-y2) /b+xl ; 

y:=(y2*a+x2) /b+yl 
end intang; 



MG 07 
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COMPUTE POINT OF INTERSECTION FROM BEARINGS 

procedure intbrg (xl ,y 1 ,bl ,x2 ,y2 ,b2 ,x,y) ; v alue xl ,yl ,bl ,x2 ,y2 ,b2 ; 
real xl,yl,bl,x2 ,y2 ,b2 ,x,y ; 

comment HUCC LIBRARY PROCEDURE MG 08: 
AUTHOR J. BOOTHROYD: 



Computes the point of intersection of the line , bearing bl, 
through (xl,yl) with the line, bearing b2 , through Cx2 ,y2) . 




This procedure uses MG 05, procedure bearing, to compute the 
bearing of the line joining (xl,yl) to (x2,y2). The base angles 
of the triangle are then computed and the point of intersection 
evaluated by the method of MG 07; 



i 



MG 09 
ifc ^ 

THREE POINT RESECTION BY ANGLES 

procedure resect (xl,yl,al,x2,y2,a2,x3,y3,a3,x,y, scale, fail) ; 

valu e xl ,y 1 , al ,x2 ,y2 , a2 ,x3 ,y 3 , a3 , s cale ; 

real xl,yl ,al,x2 ,y2 ,a2,x3,y3,a3,x,y ,scale; label fail; 

comment HUCC LIBRARY PROCEDURE MG 09: 
AUTHOR J. BOOTHROYD: 

Computes the coordinates x.'y of the point of intersection P of the 
lines P1,P2,P3 from the angles al ,a2,a3 (in radians) where al,a2,a3 
are measured clockwise from a reference observation line to the lines 
PI, P2 and P3 respectively. 




2 



The RO line may coincide with the line PI, i.e. al=0. 
The points may be in any order with respect either to themselves or 
to the resected point P. The procedure rejects all geometrically 
impossible configurations by an exit to the label fail, and special 
steps are taken to preserve accuracy in the numerically difficult 
configurations . 

The procedure is based on the formulae 

x-xl=tan bl(y-yl) 

x-x2=tan b2(y~y2) 

x-x3=tan b3(y-y3) 

b2»bl+(a2-al) , b3=bl+(a3-al) 

where bl ,b2 ,b3 are the grid bearings of the lines P1,P2,P3. 

From these may be derived the formula 

tan (b 1 ) = ( (x2-xl)cot(a2-al)-(x3-xl)cot(a3-al)+(y3"yl)) 
((y2-yl)cot(a2-al)-(y3-yl)cot(a3~al)-(x3-xl)) 

As is well known to surveyors, the case in which P,l,2,3 together form 
a cyclic quadrilateral is indeterminate. Mathematically this case yields 
tan(bl)=0/0. Near cyclic cases also produce large inaccuracies and 
are detected by evaluating the numerator and denominator of the expression 
for tan(bl) separately. The parameter scale has been provided for this 
case. The value given to scale should be twice the distance from the 
approximate centre of the knox^n points to the rnos t distant known point. 
The procedure computes 5% of this value and if both numerator and denom- 
inator are less than 0.05*scale a cyclic case is assumed; 



prQco ^uyQ rosoct<x1 ,y1 ,a1 ,x2 ? y2,a2,x3 t y3 ,a3 ,x,y , scale 9 fail) ; ffi 09 

value xl ? y1 ,a1 ,x2,y2,a2,x3 ,y3 ? a3, scale ; 

teal xl ,y1 ,a1 ,x2,y2,a2,x3,y3 ,a3 ,x,y ? scale; j.abQl fail; 

begin real t,r,eps; boolean a2zero 1 a3sero,absa2,absa3 ; switch s :=L2,L3 ,done; 
procedure solvo(xl ,y1 ,x2 ? y2,al s a2) ; 
value xl ,y1 f x2,y2,a1 ,a2; real xl ,yl ,x2,y2,a1 ,a2; 
oe^in y2:=<y2*a2~x2)/(a2~a1 ) ; 

x : =y 2* al +x1 ; y:=y 2+yl 
end sclvo; 
ep3 :=0 . r *5*scale ; 

x3:=2,3rx2; y3:=y3~y2; a3:=a3-a1 ; 
x2:=x2-x1 ; y2:=y2-y1 ; a2:=a2-al ; 
a3zer:>:=a3:=0 «€ ; a2zero:=a2=0 .0 ; 
ix a2zuro and a3zero then goto fail; 
if a3zoro then 

fie^in a2:=1 ,G/tan(a2) ; if abs(a2)> 10 5 then goto fail; 
t:=y3+y2; 

^ t ~~® »° jL^en ]begin y;=0,0; x:=a2*y2+x2-fx1 ; goto done end 
else begin al S = (x3+x2)/(y3+y2) ; goto L2 end 

^end; 

ii a2zero then 

begin a3:=1 .0/tan(a3) ; if abs(a3)> 10 5 then goto fail; 

jlf y2=0 D then begin y:=0 0; x:=a3*y3+x3+x2+x1 ; goto done ond 
else begin al :=x2/y2; goto L3 end 

end; 

a2:=l ,0/tar.(a2) ; a3:=1 «G/tan(a3) ; absa2:=abs(a2)> 10 5 ; absa3r=abs(a3)> 10 5 ; 
if alrsa2 anc absa3 jthen go to fail; 
if absa3 then 
begin t:=y3+y2; 

if t=0 # then begin y:=C»0; x : = a 2*y 2+x 2+x1 ; goto done _end 
else begin al : = (x3+x2) /( y3+y2) ; goto L2 end 

end; 

if abea2 t hen 

bugiii if y 2=0 .0 then beg in y:=0«,0; x:=a3*y3+x3+x2+x1 ; goto done end 
else beg in al :nx2/y2; goto L3 end 

ond ; 

t:=(x^*a2-<x3+x2)*a3+y3) ; r :=(y2*a2-(y3+y2) *a3~x3) ; 
if abs(t)<eps and abs(r)<eps then goto fail; 
al :=t/r; 

L2: a2: = (1 «0+a2*a1)/(a2-a1 ) ; solve(x1 ,y1 ,x2,y2,a1 ,a2) ; goto done; 
LT2 a3:=(l % 0+a3*al)/(a3-al) ; solve(x1 ,y1 ,x3+x2,y3+y2,a1 ,a3) ; 
done; end resect; 
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; MG 10 

MULTIPLE POINT RESECTION FROM ANGLES 

procedure avresect (obs , n, scale, x,y) ; value n, scale; 
real scale,x,y; integer n; array obs; 

comment HUCC LIBRARY PROCEDURE MG 10: 
AUTHOR J. B00THR0YD: 

The minimum number of known points necessary for the evaluation of 
the coordinates of a resected point is 3. Where the number of known 
points exceeds 3 some method of obtaining the best estimate of x,y 
from all computed values of x,y is required. This procedure computes 
the weighted mean of all calculations based on the selection of all 
possible combinations of 3 points selected from n points , using the 
procedure resect (MG 09) C 3 times. The array obs [ 1 :n, 1 s 4] contains 
the observational data for the n observations* For the r th observation 
X r y r \ and wei sht > should be provided respectively in obs [r ,1] , 
obs [r ,2] , obs [r,3] and obs [r,4] . For a resection based on points 
i,j,k with respective weights wti ,wtj ,wtk a group weight wti-fwtj+wtk 
is assumed. The parameter scale is used to reject near cyclic cases 
for details of which see the description of MG 09. For angular errors 
of the order of 1 second the procedure accuracy over eight (difficult) 
points is of the order of .0001 unit with a scale of 100 units. 
This means an accuracy of ,1 yd f ;x resections over distances of 
100 5 000 yds. The errors are approximately linear with angular error 
and will increase if the number of observations decreases ; 



procQdurQ inttrs(x1 ,y1 ,b1 ,x2,y2,b2 r x,y) ; v alue xl ,y1 ,b1 ,x2,y2,b2; 
real x1 ,y1 ,b1 ,x J,y2 ,b2,x,y ; 
begin real r,t; 

x2:ss3Ci2-x1 ; y2:=y2-y1 ; 

bearing(x2,y2,r,t) ; 

bl :=:1 ,0/tan(t-b1 ) ; 

t:=U+1 cO/tan(3J41592654+b2-t) ; 

x:=(x"2*b1 -y2)/t+x1 ; 

y:=r(y2*b1+x2)/t+y1 

end intbrg; 



proceduro avroooct(obs ,n, scale ,x ,y) ; valuo n, scale; 
real scale, x,y; jntogor n; array obs ; 
bogin ijit ogor i , H ,k; 

real sn,sx,ey,xl ,yl ,a1 y x2,y2,a2,xx,yy ,wti ,\vt j ,wt ; 

swr tch ?:~loop; 

sns==sx:=sy:=0 c r ; 

for gtop 1 until n do 

bg^in xl :=obGLi,1 3; yl :_obs[i ,2] ; al : =obs[i ,3] ; wti :=obs[i ,4] ; 
lor j:~i+1 step 1 until n do 

bogxn x2:=obs[j,1 ] ; y2:=obs[ j ,2] ; a2:=obs[ j ,3] ; wt j ;=obs[ j ,4] ; 
for k:=j+1 stop 1 unti l n do 

fcogin rosect(x1 ,y1 ,a1 f x2,y2,a2,obs[k,l ] ,obs[k,2] ,obs[k,3] ,xx f yy , scale, loop) ; 

wt :=wti+wt j+obs [k ,4] ; 
sn:=sn+wt ; 
sx : =sx+wt*xx ; 
sy:=sy+wt*yy ; 

loop: end k 
pi-d j 

erd i; 

ii^sn+r then begin x:=sx/sn; y *=sy/sn end 

olsg jDrint punch(3) ,££1?RESECT ERROR?, wait 

end avrosecr ; 
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THREE POINT RESECTION FROM DISTANCES 

procedure dis tresect (xl,yl ,dl ,x2 ,y2 ,d2 ,x3 ,y 3 , d3 ,x,y , f ail) ; 

v^lue Xl,yl,dl,x2,y2,d2,x3,y3,d3; 

real :;l,yl,dl,x2,y2,d2,x3,y3,d3 f x,y; label fail; 

comment HUCC LIBRARY PROCEDURE MG 11: 
AUTHOR J. BOOTHROYD: 



Computes the coordinates of a resected point P from the distances 
dl,d2,d3 from P to known points (xl,yl) , (x2,y2) ,(x3 f y3) . This 
procedure has been written to provide assistance to those making 
use of telurometery . 

2 90 
The equations are (x-xl) +(y-yl) = dl 

(x-x2) 2 +(y-y2) 2 - d2 2 

(x-x3) 2 f(y-y3) 2 - d3 2 
from which two linear equations in x and y may be obtained by 
subtraction in pairs. If the data provided results in a zero 
determinant of the two linear equations the procedure exits to the 
label fail; 
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MULTIPLE POINT RESECTION USING DISTANCES 

procedure avdis tresect (obs ,n,x,y) ; value n; 
integer n; real x,y ; array obs; 

comment HUCC LIBRARY PROCEDURE MG 12s 
AUTHOR J. BOOTHROYD : 

This procedure evaluates the coordinates x,y of a resected point 
as the weighted mean of all computations based on the use of all 
combinations of 3 points selected from n observations , using distances 
in all cases . The procedure disresect (MG 11) is called n 3 times . 
For the r th observation the data x^Jy f 4 jWeight^ should occupy 
obs [r,l] , obs [r ,2] , obs [r ,3] 5 obs [r,4] respectively of the 
array obs[l:n,l:4] . 

For a computation using observation i,j,k with respective weights 
wti,wt j ,wtk, a group weight wti+wt j+wtk is assumed; 



procjdur e distresecfc(x1 ? y1 ,d1 ,x2,y2 ? d2,x3 ,y3,d3 ,x,y,fail) ; 

value xl ,y'i ,dS ,x2,y2,d2,x3 ? y3,d3 ; 

real x1 ,y1 ,11 ? x2,y2,d2,x3,y3,d3,x P y; label fail; 

begin xeal all ? a1 2*a21 ,a22 P b*S ? b2*det ; 

al fF=(x2-xl) ; al 2:=(y2-y1 } ; 

a 21 :=(x3~x1 ) ; a22: = (y3-y1 ) ; dl s=d1 *d1 : 

bl :=<a11*(x2+x1)+a1 2*(y2+yl )-d2*d2+dl )/2 c ; 

u2:-<a^1 *(x3+x1 ) +a22*(y3+y1 )-d3*d3+d1 )/2 1 ; 

cet:=e"?l *a22-a21 *a1 2; 
if dct-=0 ,C then goto fail; 
x:=(a22*b1 -a 21 *b2) /dot ; 
y:=(a11 *b2~a21 *b1 ) /dot 
end disxresoct : 



rrocoduro aviistresoct(obs -n r x,y) ; valu e n; 
real ; intogor n; array obs ; 
b^gin int 3fe? i, j ,k; 

real sn 5 sx,sy ,x1 P y1 ,d1 ,x2 ,y2 , d2,xx r yy ,wti ,wtj ,v/t ; 

3V'i t ch s :=loop ; 

£n:=ux;=sys=( c ; 

for i:=1 step 1 until n do 

beg^n Xfl :=obn[i ? 1 ] ; yl :=obs[i,2] ; dl s=obs[i ? 3] : wti :=obs [ i ,4 ] ; 
for j:~i+1 step 1 until n do 

begir x2s=obs[ j ,1 ] ; y2:=obs[ j , 2] ; d2:=obs[ j ,3] ; wtj:=obs[ j,4] ; 
for j-fl stop 1 unti l n do 

begin distresect(x1 ,y1 ? d1 ; x2,y2,d2,obs[k 5 *l ] ,obs[k,2] ? obs[k,3] oXX,yy , loop) ; 
wt :=wti+wt j+obs[k r 4] ; 
sn:=sn+wt ; 
sx:=sx+wt*xx; 
sy:=sy+wt*yy; 

loop: end k 
end j 

end i; 

jif on*^ then begin x:=sx/sn; y 2=sy/sn end 

else print punch(3> ,££1?DISTRESECT ERROR? .wait 
end avdistresect ; 
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THE INACCESSIBLE BASE PROBLEM 

procedur e f arbase (xa,ya,xb ,yb ,aa, ab ,ba,bb ,xl ,yl,x2 ,y2) ; 

value xa,ya,xb ,yb ,aa,ab ,ba,bb ; 

real xa,ya,xb ,yb,aa, ab ,ba,bb ,xl,yl,x2 ,y2 ; 

commen t HUCC LIBRARY PROCEDURE MG 13: 
AUTHOR J. BOOTHROYD: 



Solves the problem of the inaccessible base. 




(xl,yl) 

(xl,yl) and (x2,y2) are stations whose coordinates are required, 
(xa,ya) , (xb ,yb) are stations whose coordinates are known . 
aa,ab ,ba,bb are angular observations from (xl,yl) and(x2 ,y2) on 
(xa,ya) and (xb ,yb) using the angular sign convention of procedure 
MG 07. Standard survey computational methods would compute the 
tangent of the angle of inclination of line a,b to line 1,2 and, 
from the computed bearing of 1,2 and the measured angles deduce 
the bearings of all other lines. This would, provide sufficient 
information to evaluate xl ,yl ,x2 ,y2 using procedure MG 08. This is 
not the best method if automatic computing facilities are available. 
Using the formulae described in procedure MG 07 we obtain the four 
equations :- 

xa - ( (x2-xl) cctaa-(y2-yl) ) /(cotaa + cot ba)+xl 
xb = ((x2~21) cot ab-(y2-yl))/(cot ab+cot ' " )+xl 
ya = ((y2-yl)cct aa+(x2-xl) ) / (cot aafcot ba)+yl 
yb - ((y2-yl)cot ab+(x2-2 1) ) / (cot ab+cot bb)+yl 
which, rearranged yields four linear equations in the four unknowns 
xl ,yl ,x2 ,y2 . This procedure computes the coefficients of the linear 
system and calls on LE 02 (SOLVEQ) to ob tain the solution; 



procedure f arbao e(xa ,ya ,xb *yb , aa ,ab ,ba *bb ,x1 ,y1 ,xJ?,y2> ; 

value xa , ya *xb ; yb ,aa ,ab ,ba ,bb ; 

real xa , y a ,xb , yb , aa , ab , ba ,bb ,Xl ,y1 5 x2,y2; 

begin real aaba,abbb; array A[1 s4 ,1 :5] ; 

^ " aa:=1 3/tan(aa) ; ab:=l o 0/tan(ab) ; 

bas=1 # 0/tan(ba) ; bbt-1 o 0/tan(bb) ; 

aaba;=aa+Ha ; abbbs~ab+bb ; 

A[1 ,1 ^s=At3,5]:=ba; A[2,1 ]s=A[4,2] :=bb; 

A[1 ,3]:=Ar.3,4]:=aa; A[ 2,3] :=A[4 ,4] :=ab; 

A[1 ,'Si:=^ r 2 ? S]:=A[3 f 3]:=A[4,3]:=l o ; 

Atl ,4j:=AL2 a 4]:=A[3,1 ]s=A[4,1 ]s=-1 c ; 

AD ,5] :=xa*aaba ; A[ 2 ,5] s=xb*abbb ; 

A [ b , 5 ] : ~ya*aaba ; A[4 ,5] s =yb*abbb ; 

SOL\^q(A ; 

■ 3d?=A[1 ,51; yl2=A[2,5]; x2:=A[3,5]; y2s=A[4,5] 
enc xarbaso; 
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3-DIMENSIONAL COORDINATE GEOMETRY MP VECTOR ARITHMETIC PACKAGE 

The procedures of this package are collectively designated as MG14. 
No provision is made for copying individual procedures into user's 
tapes since these are brief enough to be re-written as required. 

comment HUCC LIBRARY PROCEDURE MG14: 
AUTHOR : B. ROBINSON: 

The point P(x,y , z) in three dimensional space has one-to-one 
correspondence with the vector V which has components x,y ,z 
parallel to the three mutually perpendicular coordinate axes 
respectively, In these procedures the vector V (which may be 
visualised as the directed line-segment OP from the origin 
to the point P) is represented by an ALGOL array of three 
elements , which is declared thus 

real array V [ 1 : 3 ] ; 
The elements of this array are the components of the vector, 

V[l]«x 

V[2]=y 

V[3]«z 

In procedure calls the vector is referred to as V[j ] « 

The subscript j must be declared as an integer and the 

declaration must be valid for all blocks in which these 

procedures are employed. Inside the procedures j is assigned 

the values 1,2, or 3 as needed for operations on the components 

of the vector. Outside the procedures j should not be used . 

The procedures make extensive use of the ALGOL facility 

of calling by name. This is convenient because it makes possible 

the use of linear expressions of vectors as actual parameters in 

procedure calls, and the use of arrays of vectors or individual 

vectors. For example, a two- dimen s i on a 1 lattice of vectors with 

typical element L would be declared as 
J A mn 

real array I.[£l:ul, 12 :u2, 1:3] 



and any individual vector of this lattice could be referred to in a 
procedure call as L[m,n,j], 

The result of this is that the text of the program will be conveniently 
reminiscent of vector algebra. 



VECTOR ASSIGNMENT : 

procedure make(vi) equal to:(xj); real xj,yj; 
Examples : make(V2[j] ,Vl[jJ) ; effect is f, V2:-Vl M 

make(V3[j],Vl[j] + V2[j]); "V3:=Vl+V2 n 

make(V2[j],Vl[j]*sc) ; "V2 :=Vl*sc" 

make(M[j],(A[j]+B[j])/2.0); 
This finds the midpoint of the line-segment AB 

effect is ff M:=(A+B) /2" 

make(G[j],(A[j]+B[j]+C[j])/3.0); 
This finds the centroid of the triangle ABC, 

effect is "G:«(A+B+C) /3" 



MULTIPLE VECTOR .ASSIGNMENT : 

procedure make both(vi) and:(zj)equal to:(xj); real x j ,y j , z j ; 
Example: makeboth(V[6, j] ,temp[ j] ,V2[j ]-Vl[ j] ; 

effect is "V[6] :=temp:=V2-Vl lf 



DOT PRODUCT OF TOO VECTORS : 

real procedure dot(xj ,yj) ; real xj ,yj ; 

The procedure takes the value of the scalar product of the two vectors • 
Example: r:=dot(Vl[ j ] ,V2[ j ]) ; 

If the dot product is zero, then the two vectors are perpendicular 
or else one is a zero vector. 



LENGTH OF A VECTOR; 

real procedure length (xj) ; real xj ; 
NOTE: uses procedure dot. 

Examples : The length of the vector V from the origin is found by 

r:=length(V[j]); 
The distance between the ends of the vectors V^ and is 
found by r:«length(V2[j]-Vl[j]) ; 

COSINE OF THE ANGLE BETWEEN TOO VECTORS : 

real procedure cosang(xi «vi) escape label : (fail); 
real xj »yj ; label fail: 

The procedure takes the value of the cosine of the angle subtended 
at the origin by the two vectors. If either vector is of zero length 
then the angle is indeterminate and the procedure exits to a label, 
Examp le : cs the ta : =cos ang ( VI [ j ] , V2 [ j ] , f ai 1) ; 

To find the cosine of an angle not at the origin the call by name is 
so arranged as to simulate a change of origin to the apex of the angle* 
Thus where the angle /ABC is represented by the three vectors a[j], 
b[j] and c[j], 

cstheta:«cosang(a[j]-b[j], c[j]-b[j], fail); 
NOTE: uses procedures dot and length, 

VECTOR OF UNIT LENGTH PARALLEL TO A GIVEN VECTOR : 
procedure unit(nj)parallel to :(xj) escape label: (fail); 
real nj,xj; label fail; 

Example: unit(V[ j] ,V[ j3 ,fail) ; converts the vector V into 

a unit vector in the same direction. The procedure will exit by the 

label if the vector given to it has zero length. 

NOTE: uses procedures dot and length 



CROSS PRODUCT OF TOO VECTORS ; 

procedure cross(xj,yj) is written into;(zj); real xj ,yj , zj ; 
Example; cross(Vl[j], V2[j], V3[j]); effect is "V3;=VlxV2" 
The direction of V3 is perpendicular to the plane in which VI and V2 lie. 
If VI and V2 are parallel vectors then V3 will be a zero vector. If VI 
lies along the x-axis and V2 lies along the y-axis, then V3 will lie 
along the z-axis in accordance with the convention, this may be expressed 
by the statement that right-handed coordinates are used (remember that 
VlxV2*-V2xVl) . 

VECTOR OF UNIT LENGTH PERPENDICULAR TO TWO GIVEN VECTORS ; 

procedure perpunit(ni) perpendicular to;(xj,yj) escape label; (fail) ; 

real nj,xj,yj; label fail; 

Example ; perpunit (n [ j ] , VI [ j ] , V2 [ j ] ) ; 

NOTE; uses procedures dot, length, unit and cross. 

BOX PRODUCT OF THREE VECTORS ; 

real procedure box(xj ,yj ,zj) ; real xj,yj,zj; 

The procedure takes the value of the scalar triple product of the 
three vectors. If the three vectors are coplanar with the origin 
then the product is zero. It can therefore be used to test whether 
any four points are coplanar, by evaluation of 

q:*box(a[j]-b[j], c[J]-b[j], d[jj-b[j]); 
NOTE; uses procedures dot and cross. 
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comment vector procedures; 



procedure make(yj,xj> ; real xj,yj; for j:=1,2,3 do yj:=xj; 

procedure makeboth(yj,zj,xj> ; real xj,yj,zj;' for j:=1,2,3 do yj:=zj:=xj; 

real procedure dot(xj,yj); real xj,yj; 

begin real t; t:=0.0; for 3: =1 ,2, 3 do t :=t+xj*yj; dot:=t 
end dot product; 

real procedure length(xj) ; real xj; length: =sqrt(dot(xj,xj)) ; 

real procedure cosang(xj,yj, fail); real xj.yj; label fail; 

begin real d; d:=length(xj)*length<yj) ; if d=0.0 then goto fail; 

cosang:= dot(xj,yj)/d 
end cos angle; 

procedure unit(ni.xi.fail) ; real nj.x.i; label fail; 

begin real d; d:=length(xj) ; if drO.O then goto fail; 

for j:=1,2,3 do nj:=xj/d 
end unit; 

procedure cross(x.i.vi.zi) : real xj,yj,zj; 
begin real x1 ,x2,x3,y1 ,y2,y3; 
3 : =1 ,* xl : =xj ; yl : =yj ; 
j:=2; x2:=xj; y2:=yj; 
j:=3; x3:=xj; y3:=yj; 

Zj:=x1*y2-x2*y1; j:=2; zj : =x3*y1 -xl *y3; 
3 ' =1 ,* z 3 : =x2*y3-x3*y2 
end cross product; 

procedure perpunit(n.i.x.i.vi.fail) ! real nj t xj,y.j; label fail; 

begin array a[1:33; cross (xj,yj, at j]> ; unit(nj,a[ j],fail) end ; 

real procedure box(xj,yj,zj) ; real xj,yj,zj; 

begin array a[1 :33; cross <xj,yj,a[ j]) ; box:= dot(a[j],zj) 
end box product; 
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COMPUTE FOURIER COEFFICIENTS 

procedure fourier(f ,a,b ,n) ; value n; integer n; real array f ,a,b; 
comment HUCC LIBRARY PROCEDURE MHOl: 
AUTHOR: J. BOOTHROYD 

Evaluates the coefficients of the Fourier expansion of 
a function f(x) specified at the n sample points 
f[0],f[l], • f [n-l]« n nay be odd or even and 

on exit from the procedure the arrays a,b[0:n div 2] 
will contain the coefficients appropriate to either case 
as follows 
n even (n ~ 2N) 

N-1 

f(x)-a[0] + S(a[k]*cos(kirx/N)+b[k]*sin[kTTx/N))+a[N]cosirx 
k=l 

with b[0] = b[N] - 0. 
n odd (n =» 2N+1) 

N 

f(x)=a[0] + E(a[k]*cos(27rkx/n)+b[k]*sin(2T:kx/n)> f 
k=l 

with b[0] - 0. 
The algorithm uses the recurrence relations described 
in R« W« Hamming - Numerical Me thods for Scientists and 
Engineers* page 72, 
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procedure fourier(f ,a,b,n) ; value n; integer n; real array f ,a,b; 
begin integer ndiv2,k,m; real c,s f vk,v1 ,u0,u1 ,um,t ,ck,nby2; 

ndiv2:= n div 2; nby2:= n/2; t:= 6.2831853/n; 

c:= cos(t) ; s:= sin(t) ; vk:= 0.0; vl : = -1 ; 



begin t:= c*vk; ck:= t-vl; vl := vk; vk:= ck+t; 
t:= ck+ck ; ui:=0.0; um: = f[n-l]; 
for m:= n-2 step -1 until 1 do 

begin u0:= ut ; ul := um; um:= t*u1-u0+f [m] end ; 
a[k]:= (ck*um-u1 +f [0] )/nby2 ; 
b[k] := s*vl*um/nby2 



end ; 

a[;]:= a[ 3/2.0; if ndiv2*2=n then 
begin afndiv23:= a[ndiv2]/2, 0; 
b[ndiv23:=0.0 

end 

end fourier; 



for k:= step 1 until ndiv2 do 




» 



MHO 2 



.1 



SUM FOURIER SERIES 



real procedure sunf ourier(a,b,n,x) ; value n,x; real x; 
integer n; real array a,b; 

consent HUCC LIBRARY PROCEDURE MH02: 
AUTHOR: J, BOOTHROYD 

Evaluates the function f(x) for given argument x 
by sunning the Fourier series of the function, the 
coefficients of which occupy arrays a,b[0:n div 2] 
as a result of using procedure MH01 t 



real procedure sumfourier(a, b,n,x) ; value n,x; real x; integer n; real array a # b; 



comment evaluates the function f (x) for given argument x by summing the Fourier series whose 
coefficients are the elements of a,b[0:ndiv2] derived by the use of procedure fourier, i.e. 
n even (n=2N) 

f(x)= a[0]+sigma(a[k3*cos(pi*k*x/N)+b[k]*sin(pi*k*x/N>,k,1 f N-l) + a[N]*cos<pi*x) b[0]=b[N]=0 
n odd <n=2N+1) 

f(x)=a[3]+sigma(a[k]*cos(2*pi*k*x/n)+b[k]*sin(2*pi*k*x/n) f k # l # N) , b[0] = 0; 
begin real vk,v1 ,sum,c f s,ck; integer k,ndiv2; 

ndiv2:=n div 2; x:=6. 283l853*x/n; c:=coe(x); s:=sin(x); 
sum:=a[o]; vk:=l # 0; vi:=o,0; 
for k:=l step 1 until ndiv2 do 
begin x:=c*vk; ck:=x-v1 ; 
vl : =vk ; vk : =ck+x ; 
sum : =sum+a C k ] *ck+b [k]*s*vl 

end ; 

sumf ourier : =sum 
end sumf ourier; 
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LOCATE ROOT OF f (x)=Q 

real procedure bisec(f ,x»a,h,eps , error) ; value a,b,eps ; 
real f ,x,a,b,eps ; label error; 

comment HUCC LIBRARY PROCEDURE MROl: 
AUTHOR: NOT KNOWN 

A procedure to find a root cf f (x)=0 in the interval 
a to b by continued bisection. If f (a) and f(b) 
have the same sign the procedure exits to the label 
error. The error escape facility may be used to 
determine the interval (a,b) wi thin which some root 
of interest lies. Suppose that f (x) has several 
positive roots X} x 2 X3 and it is known that the 

root seperation (x 2 - xi) , (X3 - x 2 ) etc. is at 
least z but the location of xi is in some doubt f 
The routine 

xa := -z; 
strobe; xa := xa + z: 

y := bisec(f (x) ,x,xa,xa+z ,10" 4 S strobe) 
will increment xa until f (xa) and f (xa+z) have 
opposite signs when the bisec routine will evaluate 
the required root to an accuracy of 10 




j MR01 



real procedure bisec(f,x, a,b,eps,error) ; value a,b,eps; real f,x,a,b,eps; label error; 
comment the procedure finds a root of f (x)=0 in the interval a to b by- 
continued bisection. If f (a) ancf f (b) have the same sign the procedure exits 
through the label error; 

begin real q; x:=a; q:=f; x:=b; 
if f*q>0 then goto error; 

q : =(b-a) *si gn( q) /2 . ; x : =( a+b) /2.0; eps : =eps/2 . ; 
£or q:=q/2.0 while abs(q»eps do xj=x+(if f>0 then q else -q) ; 
bisec:=ac 
end bisec; 



) 
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GENERAL SUM SERIES 

real procedure sigma( tk,k,a,b) ; value a,b; integer k,a,b; real tk; 

comment HUCC LIBRARY PROCEDURE MS01: 
AUTHOR: J. EOOTHROYD 

A Jensen procedure to evaluate the sum over k from 
a to b of tk with positive unit increments in k. 
The parameters tk and k are called by name. As examples 
of the flexibility of this procedure the calls, 

1. sigma(a+i*d t i^0 > 3) 

2. sigma(a*rti,i>0,n-l) 

3. sigma(a[i}*b[i],i,l,N) 

will evaluate respectively 

1. The sum of 4 terms of an arithmetic progression 

with first term a and common difference d, 
2 # The sum of n terms of a geometric progression. 
3, The inner product of the vectors a,b[l:N] . 

These examples illustrate the dependence of the actual 
first parameter on the variable actually used as the second 
parameter* 



MS01 



real procedure sigma(tk,k,a,b) ; value a,b; real tk; intege r k,a,b; 
comment the sum over k from a to b of tk, with positive unit increments in k; 
begin real sum; 
sum : =0 . ; 

for k:=a step 1 until b do sum:=sum+tk; 
sigma: =sum 
end sigma; 



t 
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20 ROMBERG INTEGRATION 



real procedure rombergl (x, fx, a,b,order ,eps) ; value a,b,eps,order; 

real x,fx,a,b,eps ; integer order; 
begin integer i f j,mark; real ti,k,h,n,m,t,fac,del, abs t m , abs t o 1 , abs t i ; 

array aa[1 ; order] ; switch s : =test , newi , done ; 

x:=a; ti:rfx; x:=b; ti:=aa[1 ] :=(ti+fx)/2 # 0; 

k:=h:=b-a; 

if k=G # Q then begin i:=0; goto done end ; 
n:=1 # 0; abstol:= 10 -8/abs(k) ; 
for i : step 1 until order do 
begin h:=di/2.0; m:=0.0; 

for x : =a+h step k until b do :=m+fx; 

m:=m/n; t:=ti; j:=i; 

absti :=abs(ti) ; 
test: abstm: =abs(m-t) ; 

if ti+O.O then 

begin if abstm/absti<eps then goto done end 
else 

begin if abstnKabstol then goto done end ; 

if i j then 

begin if mark=0 then 

begin t : =aa[ j ] : =m+(m-t ) /f ac ; 

mark : =1 ; del:=3.G*(fac+l .0); 

fac:=fac+del 

end 
else 

begin j:=j-1 ; if j=0 then goto newi; 
m: =t ; t:=aa[j]; mark:=0 

end ; 

end 
else 

100 begin mark: =0; fac:=3 # 0; del:=0 # 0; 

ti:=m:=aa[ j] :=(m+t)/2.0; j:=j-1 

end; 

goto test ; 
newi : k:=h ; n : =n+n 
end i ; 
i:=0; 

done: r omberg : =( b -a ) * ( i f i>0 then (t+m) /2.0 else aa[ 1 ] ) 
end romberg; 
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ADAPTIVE SIMPSON INTEGRATION 
real procedure sinps (f , x, a, b, eps); value a, b, eps; real 
f, a, b, eps; 

comment HUCC LIBRARY PROCEDURE MS02 : 
AUTHOR J. BOOTHROYD : 

A modification of Algorithm 233 COMM. A. CM. VOL 7 NO. 6 
JUNE 1964 p 348. The procedure evaluates the integral of 
f with respect to x from lower bound a to upper bound b. 

b 

simps : * / f dx 
a 

The parameter eps determines the acceptable relative error between 

successive evaluations of the integral I. If 1^ and I ^ are 

successive approximations, the process will terminate when 

abs ((r - 1^ > / y < eps 

unless I = in which case, to avoid division by zero, the 

process terminates when the absolute difference between successive 

values of I differ by less than 10 ~ 8 . As an example of the use 

of simps, the call 

y : - simps ( i.o/x, x, 1.0, B, .0001) 

B 

will assign to y the value r , ly v - , r „, 

{(Vx) dx = in [B] 

with a tolerance of .0001. 

The procedure operates correctly for the cases a < b, a ■ b and a > b. 



real procedure simps(f ,x ? a, b, cps) ; value a,b, eps ; r^eal f ,x,a,b,, eps ; 
begin real Z1 , Z2, Z3,h p k; switch s := again, done ; 

if a=b then b egin h:=0. 0; poto done end; 

x:=a; Zl:=f; x:=b; Z1 :~Z1+f ; k:=(b«a)/2*Q; 

x:=a+k; 22: =f; Z3:=Z1-:-4,,0*Z2: Z1 :=Z1+2 G*Z2j 
again 2 Z2:=0; h:=k/2; 

for x:=a+h ste p k until b do Z2 : =Z2+f ; 

Z1:=Z1+4.0*Z2; 

if ZlfcO then b egin if abs( (Z1-2. 0*Z3)/Z1 )<eps then goto done end 

els e if abs(Z1-2.0*Z3)*h/3 9 0<i -3 then goto done ; 
Z3:=Z1; Z1:=Z1-2.0*Z2; k:=h; goto again; 
done : simps : =h*Z1 /3. 

end simps; 



1 
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EVALUATE DEFINITE INTE GRAL 
real procedure havie (x, a, b,eps, integrand ,m,ma$k) ; 
, v , a l u £, a,b,eps,m,mask; integer m; real a ,b , eps , integrand ,x 5 rna>k ; 
comment HUCC LIBRARY PROCEDURE MS 03 
Author : ACM 257 

Performs numerical integration of integrand with respect 
to x over the interval a to b 

b 

havie :« y integrand dx 

eps is the convergence criterion, m is the maximum order 

\ 

of approximation to be considered in attempting to satisfy 
the eps criterion and mask is a value supplied by the user 
which will be returned by havie If convergence is not obtained, 
The method is basically the trapezium formula with higher 
order corrections (m specifies an upper limit to the number 
of these performed) „ 



Examples of calls are :~ 

ir/2 

(a) To evaluate y- j cos x dx 



y:= havie (x ? 0.C 3 l e 5707963 ? 10" 6 ? ccs(x) ,10,9.99999) ; 



(b) To evaluate z= ; £~ y2 dy 



4.: 

/ 





z:= havie(y,o.O,4.3,10 6 ,exp(y*y) , 10 ,9 . 99999) 



real procedure havie(x,a f b,eps, integrand, m, mask) ; value a $ h 9 eps,m f mask; integer in; real a, b, x , eps , integrand, mask 
begin real h , endp 1 3 , sumt f sumu 9 d ; i nteg er i, j,k,n; swit ch s:=estimate, test, exit ; 
real array t,u f tprev,uprev[ 1 :m] ; 

x:=a; endpts : -integrand ; x J =b ; endpte:=0,5*(integrand+endpts) ; 
sumt :=0.0; i:=n:=1 ; h: =b-a; 
estimate: t[ 1] :=h*(endpts+sumt) ; sumu:=b,0; x:=a-h/2«0; 

for j:=1 step 1 unti l n do begin x:=x+h; sumu :=sume+integr and end ; 
uC 1 ] : =h*sumu ; k : =1 ; 

test: if abs<t[k]-u£k]> < eps then begin havie:=0.5*(t[k]+u[k]) ; goto exit end; 

if k = i then 

begin d: =2. Of <k+k) ; t [k-f i ] : =(d*t[k] -tprev[k])/(d~l c 0) ; tprevCk] :=t[k] ; 
u[k+ 1 ] : =(d*u[k] -uprev[k] )/(d-1 „ 0) ; uprevfk] :=u[k] ; k:=k+1 ; 
' _if k=ia then begin havie:=mask; goto exit ond; goto test 

end; 

h: =h/2 # 0; sumt :=sumt+sumu; tprevCk] :=t[k] ; uprev[k] :=u[k] ; 
it=i+1 ; n:=n+n; goto estimate; 
exit: end havie; x 



> 
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LAGRANGE INTERPOLATION 

real procedure Lagrange (x,y,arg,n,m) ; value arg,m,n; 
array x,y ; real arg; integer m,n; 

comment HUCC LIBRARY PROCEDURE MT01: 

AUTHORS: J. BOOTHROYD, G. WALKER. 

The array y [0:n] contains function values y(x) 
at the sample points of x[0:n] which is sorted 
in ascending order. The procedure finds an 
approximation to the function at the specified 
point arg by evaluating an mth order polynomial 
(m<n) , choosing the appropriate subset of m+1 
sample points so that these are evenly distributed 
about arg. If the requested m exceeds n 
the assignment m:=n occurs . 

For arg < x[0] the procedure extrapolates using 

x[0] x[m] 

For arg > x[n] the procedure extrapolates using 

x[n-m] .... x[n] 
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real procedure Lagrange(x,y,arg f n,m) ; 
value arg,m,n; array x,y; real arg; integer n,m; 
begin integer i, j,min; 

real term,f ac, yest ; switch SS:=out; 
integer procedure setmin(L) ; label L; 
begin integer i ; 

switch S:=found; 

for i:=0 step 1 until n do 

if arg<x[i] then goto found; 

i:=n; 

found: j^f argrac[i] then begin yest:=y[i] ; goto L end ; 

i:=i-m div 2-1 ; 

setmin: = if i<0 then 3 else if i+m>n then n-ra else i 
end setroin; 
if m>n then m:=n; 
min:= setmin(out) ; 
fac:=1.G; 

for j:=min+m step -1 until min do 

fac:=fac*(arg~x[j]); 

yest:=:0; 

for i:=min+m step -1 until min do 

begin term:=y[i]*f ac/(arg-x[i]) ; 

for j:=min4m step -1 until i+1,i-1 step -1 until min do 
term:=term/(x[i]-x[ j]) ; 
yest ; =yest+term 

end ; 

out : Lagrange : =yest 

end Lagrange; 



i 
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AITKEN UNEQUAL INTERVAL INTERPOLATION 

Ok 

real procedure ^titken (x > y,arg s n,m ,s ; value arg^rum; 
array x,y ; real arg; integer m, n; 

comment HUCC LIBRARY PROCEDURE MT 02; 
Author : J, Boothroyd. 

This procedure performs the same function as MT 01 

but is more efficient and economical. The array y[0:n] 

contains function values y(x) corresponding to the sample 

points of x[0:n] which is assumed sorted in ascending order. 

The procedure estimates the value of the function corres- 

ponding to a specified point arg by evaluating an m LU order 

polynomial (m<n) , choosing the appropriate subset of nrfl 

sample points evenly distributed about the value of arg. 

If the requested value m exceeds n the assignment m:~n occurs . 

For arg < x[0] the procedure extrapolates using x[0] x[m] . 

For arg > x[n] the procedure extrapolates using x[n-m] . . .x[n] ; 



MT02 



real procedure aitken<x,y,arg,n,nO ; value arg,m,n; SSESL **'. i2£SS22 

ffiirT integer i^nflessi; real fi.ai; real arrav. z,f[0:ml; switch .tseut. 

integer procedure setniin(L) ; labe l L; 

begin integer i; switch s:=fouhd; 

' — fSTir^O step l until n do if arg<x[i] then goto found; 

x * ~n * 

found: if arg=x[i] then begin f[u]:=y[i); goto L end; 

i:=d-m div 2-1 ; 

setmin:= If i<0 then else if i+n>n then n-n else i 
end setmin; 

if m>n then mi =n; j: =setmin(out) ; end , 
fSr i:=07tep 1 until m do 'begin z[i] :=arg-x[ j] ; f[i]?=y[jl; Jt«J+1 eud, 

mlessl:=m-1 ; 

fcr i;= step 1 until mlessl do 
begin fis=f [il; si:=z[i]; 

for j:= i+1 step 1 until m do 

f[j]:= fi+zi*(fTj]-fi)/(2i-2Ul) 

end; 

out: aitken:=f [m] 

end ; 



l 
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aItken equal interval interpolation 

i£3l procedure equipol (xbase, y ,arg,n f m f h) ; value xbase ,arg,m,n,h; 
real xbase, arg, h; array y; integer m,n; 

comment HUCC LIBRARY PROCEDURE MT 03 
Author : J. Boothroyd. 

An equal interval version of MT 02, Array y[0;n] contains 
values of a function y(x) corresponding to the equal 
interval values of the argument xbase, xbase + h, xbase + 2h. . . 
«...., xbase + nh. The procedure estimates the value of the 
function for some specified value arg by evaluating an m th 
order polynomial (m£n) 5 choosing the appropriate subset of 
nri-1 sample points evenly distributed about the value arg. 
If m exceeds n the assignment m:=n occurs . 

For arg < xbase the procedure extrapolates using y [0] . • . .y [m] . 
For arg > xbase+nh the procedure extrapolates using 
ytn-m] y[n] ; 



MT 04 

AITKEN N th Q VT^v T T ^ zcr > m TCW ~ 

real procedure ait(z,f ,arg,n) ; value arg,n; integer n; 
real arg ; array z,f ; 

comment HUCC LIBRARY PROCEDURE MT 04 
Author : J. Boothroyd 

An abridged version of MT 02. Array f [0:n] contains values 
of a function corresponding to the values of the argument in 
z[0:n] which may be in any order. The procedure estimates 
the value of f (z) corresponding to some specified value of 
z=arg by evaluating an n th order polynomial, using all the 
values of z,f [0:n] ; 
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real procedure equipol<xbase,y,arg,n,m,h) ; value xbase , arg , m , n , h : real xbase, arg,h; array y; integer m, 
begin integer i^.mlessi; real jh f fi; array f[0;m]; — 
if m>n then m:=n; i:=entier((arg-xbase)/h)-m div 2; 
j:= if i<0 then else if i+m>n then n-m else i; 
for i;= step 1 until m do f [i] : =y[i+j] ; 
<^g£=arg^^&5 
xnlessi : =an-1 ; 

for i:=0 step 1 until ralessl do 
begin fi:=f[i]; jh:=h; 

for j:=i+1 step 1 until m do 

begin f[j]:=fi+arg*(f[j]-fi)/jh; jh:=jh+h end ; 
arg:=arg-h 

end; 
equipol :=f [m] 
end equipol; 

MT04 

real procedure ai t(z , f f arg , n) ; value arg,n; integer n; real arg; array z,f; 

begin integer i, Xnlessi ; real fi,zi,u; nlessl : =n-1 ; 
for i:=0 step 1 until nlessl do 

begin fi:=f[i]; zi:=z[i]; u:=arg-zi; 
for j:=i+1 step i until n do 
f [ j] : =f i+u*( f [ j] ~f i) /(z[ j] -zi) 

end; 
ait:r:f[n] 

end; 



SECOND ISSUE 
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ESTIMATION OF DERIVATIVE OF GIVEN FUNCTION - UNEQUAL INTERVALS 

real procedure dydx(x,y ,arg,n,m,es t) ; value arg,m,n; 
array x,y ; real est , arg; integer m,n; 

comment HUCC LIBRARY PROCEDURE MT 05 
Author : J. Boothroyd. 

Array y [0 :n] contains values cf a function y(x) corres- 
ponding to the sample values of the argument in x[0 :n] which 
is assumed sorted in ascending order. The procedure estimates 
the value of the derivative of y (x) at some specified value 
of x=arg by evaluating the derivative of an w~ order 
polynomial (m<n) 9 choosing the appropriate subset of mfl 
sample points evenly distributed about the value of arg. 
If m>n the assignment m:~n occurs. 

For arg <x[0] the procedure extrapolates using x[0] . . • .x[m] 
For arg >x[m] the procedure extrapolates using x[n-m] • , . .x[n] . 
The output parameter est delivers an estimation of y(arg) 
yielding a value which is the same as would be obtained 

from the use of MT 02. 

In so far as derivative estimation is a numerical process 
susceptible to large errors this procedure should be used 
with caution using values of m not exceeding 5; 
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real procedure dydx(x,y ,arg,n,m,eot) ; value arg, in, n; array x,y; real arg,est ; integer m,n; 
bggin integer i,j, mlessl; real f i , zi ,dif f i , zizj, f jf i ; real array z, f ,dif f [G:m] ; 

integer procedure setmin; 

begin integer i ; switc h s:=found; 

for i:=0 stop 1 until n do if arg<x[i] then goto found; 
i :=n; 

found: i : =i -m div 2-1; 

setmin: = i£ iO then Q else if i+m>n then n-m else i 
end setmin; 
if m>n then m:=n; j:=setmin; 

for i:=0 step 1 until m do begin z[i] :=arg~x[ j] ; f[i]:=y[jl; dif f [i] :=0.0; j:=j+1 end; 



mlessl : =m-i ; 

for i:= G step 1 until mlessl do 
begin fi:=f [i]; zi:=z[i] ; dif f i :=dif f [i ] ; 
for j:= i+1 step 1 until m do 
begin zizj:=zi-z[ j] ; fjfi:= f[j]-fi; 

diff[ j]:=diffi+(f jfi+zi*(diff [ jl-diffi))/zizj; 
f[j]:= fi+zi*f jfi/zizj 

end 

end; 

est :=f [m] ; 
dydx:=diff [m] 

end dydx; 
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ESTIMATION OF DERIVATIVE OF GIVEN FUNCTION EQUAL INTERVALS 

real procedure equidydx (xbase ,y,arg,n,ni, h, est) ; value xbase, arg ,m,n,h; 

real xbase, arg, est, h; array y; integer m,n; 

comment KUCC LIBRARY PROCEDURE MT 06 
Author : J. Boothroyd, 

An equal interval version of MT 05. Array y [0:n] contains 
values of a function y(x) corresponding to the sample values 

of the argument xbase, xbase+h, ,xbase+nh. The procedure 

estimates the value of the derivative of the function for some 

+*Vi 

w specified value x=arg by evaluating the derivative of an m 

order polynomial (m<n) choosing the appropriate subset of nri-1 
sample points evenly distributed about arg. If m>n the 
assignment m:=n occurs. 

If arg < xbase the procedure extrapolates using y[0]....y[m] 

If arg > xbase+nh the procedure extrapolates using y [n-m] . .y[n] . 

The output parameter est delivers an estimation of the function 
value y(arg) and yields a value which is the same as would 
be obtained from MT 03. 

This procedure should be used wi th discretion and for values 
of m not exceeding 5; 



( c 
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re* 1 procedure equidydx(xbase,y,arg,n,m,h,est) ; value xbase,arg,m,n,h; real xbase,arg,est,h; array y; integer m n- 

begin integer i, j,mlessl ; real jh,fi,diffi,f jfi ; array f,diff[0:m3: ' ' 

if m>n then m:=n; i:=entier((arg-xbase)/h) -m div 2; 
j:= jlf i<Q then else if i+m>h then n-m else i; 

for i:= step 1 until m do begin f [i] :=y[i+j] ; dif f [i] :=0 # end ; 

arg:=arg-rj*h; 

mlessl : =m-1 ; 

for i:=0 step 1 until mlessl do 

begin fi:=f[i] ; jh:=h; dif f i : =dif f [i] ; 
for j:=i+1 step 1 until m do 
begin f jfi:=f[ j]-fi; 

diff[j]:=diffi+(f jfi+arg*(diff[j]-diffi))/jh; 
f[j]:=fi+arg*f jfi/jh; jh:=jh+h 

end ; 
arg:=arg-h 

end; 
est:=f[m]; 
equidydx: =dif f [m] 
end equidydx; 



SECOND ISSUE 
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REARRANGE ELEMENTS OF VECTOR 



procedure ?ERMB(b,r,n) ; value n; real array b; integer array r; 
Integer n; 

comment HUCC LIBRARY PROCEDURE NC01 : 
AUTHOR J, BOOTHROYD : 



A procedure which rearranges the elements of b[l:n] so that 
b[i] := b[r[i]] i - l,2, Bi , t n« The procedure should be used 
after solving Ax = b by the use of LEO 3 and LE04 in those cas 
where further processing of the solutions x^ requires x i to 
occupy b[i] i = l y 2,»« l# , «n; 



proc edure PERMB(b p r,n) • valua n; roni array b; integer a rray r; jj^£S p £ 

comment rearranges tho elements of b[1:n] so that b[i 1 j:-b[r[i] 1 r inl , 2? c a 

b egin integer i , ■ k ; £eal w; switch S:= L; 
for n step -1 until 2 do 
begin k: = r[i] ; 

L: if k£i then 

begin if k>i then beg in ks=r[k] ; goto L end; 

w:= b[i]; b[i]:= b[k] ; b[k3:= « 

end 

end 

end PERMB; 



? 
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PERMUTE ROWS OR COLUMNS OF MATRIX 

comment mxperm(a,b, j ,k,s »d,n,p) ; value n; real a,b; integer array s,d; 
integer j,k,n,p; 

comment HUCC LIBRARY PROCEDURE NC02: 
AUTHOR: J, BOOTHROYD 

A procedure using Jensen 1 s device which exchanges rows 
or columns of a matrix to achieve a rearrangement 
specified by the permutation vectors s ,d[l:n] . Elements 
of s specify the original source locations while 
elements of d specify the destination locations. 
Normally a and b will be called as subscripted variables 
of the same array. The parameters j ,k nominate the 
subscripts of the dimension affected by the permutation, 
p is Jensen 1 s parameter. As an example of the use of 
this procedure suppose r,c[l:n] to contain the row and 
column subscripts of the successive matrix pivots following 
an in-situ matrix inversion. The two calls 

nxperm(a[j ,p] ,a[k f p] ,j »k,r,c f n,p) 
and mxperm(a [p ,j] ,a[p 9 k] , j ,k,c f r,n f p) 

will respectively perform the required rearrangement of 
rows and columns. 
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procedure mxperm(a,b, j,k,s,d,n,p) ; value n; real a,b; integer array s,d; integer j,k,n,p; 
begin integer array tag,loc[l:n] ; integer i,t,tagj; real w; 

comment set up initial vector tag number and address arrays; 
for i:=1 step 1 until n do tag[i] :=loc[i] :=i ; 
comment start permutation; 
for i:=T step 1 until n do 

begin t:=s[i] ; j:=loc[t]; k:=d[i]; 
if j* k then 

begin for p:=i step i until n do 

begin w:=a; a:=b; b:=w end ; 
tag[ j] :-tag[k] ; tagfk] :=t ; tagj:=tag[ j] ; 
loc[t] :=loc[tagj] ; locftagj] :=j 
end jk conditional 

end i loop 

end mxperm; 
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PERMUTATION OF ELEMENTS OF A VECTOR 

procedure vectorperm(m,d,n,mode,endperm) ; value n,mode; integer n,mode 
array m; integer array d; 

comment HUCC LIBRARY PROCEDURE NC03: 
AUTHOR: J. BOOTHROYD 

A procedure for generating, at each entry, a permutation 
of the elements m[l] ,m[2] . . # . m[n] of m[l:n] . (n-1) ! 
successive entries generate all n! permutations. 
The permutation is controlled by a variable radix 
counter, the array d[2:n], with digit positions d[2],d[3] ... d[n] 
where the subscript indicates the radix value. Starting with 
d=(0,0,0 .... 0) one is added to the counter at each entry to 
the procedure. One and only one element of d increases each 
time and all element positions below this are reset to zero. 
Denoting by k the subscript position of the counter element 
which increases the permutation rules are:- 

(k odd) or (k even and d[k]<2) interchange m[k],m[k-l] 
k even for 2<d[k]<k interchange m[k] ,m[k-d[k]] 

A call of vectorperm initialises d to (0,0,0,0, .. .0) with mode=l 
preparatory to further calls with mode=2. At the n! th call the 
counter d resets to zero and the procedure exits to the label 
endperm; 



procedure vectorperm(m,d f n f mode,endperm) ; value n,mode; integer n,mode; 
integer array d; array m; label endperm; 

begin integer k, j; real temp; switch s : =set , run, swap , exit ; 
goto s[mode] ; 

set: for k:=2 step 1 until n do d[k]:=0; goto exit; 

run: j:=-1; 

for k:=2 step 1 until n do 

begin if d[k]:{Jc~l then goto swap; d[k]:=0; j:=-j end; 

goto endperm; 
swap: d£k]:=d[k]+1; 

if j*1 then j:= if 2>d[k] then 1 else d[k] ; 

temp:=m[k] ; m[k] :=m[k-j] ; m[k-j]:=temp; 
exit: end vectorperm; 



NC03 



NC04 



PERMUTE ROWS OR COLUMNS OF MATRIX 

procedure mxperml(a,b, j ,k,s,d,n,p) ; value n; real a,b; integer array s,d; 
integer j,k,n,p; 

comment HUCC LIBRARY PROCEDURE NC04: 
AUTHOR: J, BO0THR0YD 



A procedure which is the same as, but avoids the use of the 

local arrays of, NC02. It is, therefore, more economical of 

space but slower in operation. The Identifiers of NC02 and 

NC04 agree over the first six characters, the parameter specifications 

are identical and NC04 is thus a direct replacement for NC02. 



NC04 



procedure mxpenn1(a,b f j,k f s,d,n,p) ; value n; real a,b; integer array s,d; integer j,k f n,p; 
begin integer i; switch ss:= scan; real w; 
for i:= n step -1 until 2 do 
begin j:= s[±]; k:= d[i] ; 
scan: i/£ jik then 

begin for p:= i+1 step 1 until n do 

if j=d[p] then begin j:= s[p]; goto scan end ; 

for p: = 1 step 1 until n do begin w: = a; a: - b; b:= w end 

end 

end 

endmxperm ; 



NC 05 



PRE OR POST MULTIPLY MATRIX BY PERMUTED IDENTITY MATRIX 
procedure permx (a,b , j ,k,r ,n 3 p,inv) ; value n,inv; 
real a,b; integer j ,k,n,p ? inv; integer array r; 

comment HUCC LIBRARY PROCEDURE NC 05 
Author : J. Boothroyd 

A procedure using Jensen's device whereby a matrix A[l:n,l:n] 
may be either pre or post multiplied by either a permuted 
identity matrix Ir or by the inverse of Ir. All row (or column) 
exchanges are performed on A without the use of an array of the 
same size . The inverse of Ir is its transpose and both Ir and 
* r ^ may be defined by a permutation vector r [1 :n] as follows 

i*j 1 i*j 
Ir[r[i],j] - Ir [i,r[j]] = 

1 1-J 1 i=j 

for i=l,2, n, j=l,2, ,n. 

Values of the parameter inv accomplish multiplication by Ir. 
Otherwise, for inv < (inv = -1 is appropriate) multiplication 
is by Ir . 

The following calls illustrate how to achieve pre or post 
multiplication 

Premultiplication by Ir permx(A[ j ,p] ,A[k,p] , j ,k,r,n,p,0) 
Post multiplication by Ir ^ perrax(A[p , j ] ,A[p,k] , j ,k,r,n,p ,-1) ; 



( < 
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procedure permx(a 9 b, j,k , r ,n,p, inv) ; value n,inv; reg.1 a,b; intege r j 9 Ic,n 9 p,inv; ijj^gg®£ 21 r J 
begin integer i ; real w ; 
procedure itori ; 
begin switch s:=L; ' 

for i:=n £top -1 until 2 do 
begin k:cr[i] ; j : =i ; 

b egin if k> j then be.^iii k:=r[k] ; goto L end end; 
for p:=1 step 1 until n do 

begin v;:=a; ~a:=b;"'~"b:=:w end 

end 

end itori; 
procedure ritci ; 
begin switch s:=sean; 

for i : ~n step -1 until 2 do 
begi n j:=i ; k:=rfi]; 
scan; if jik then 

begin for p:=d+1 step 1 until n do 

if j=r[p] then begin j:=p; goto scan end; 
for p:=l step 1 unti l n do 

begin v;:=a; a;=b; b:=v; end 

end 

end 

end ritoi ; 

if inv>0 then itori els e ritoi 
end permx; 



NC 06 j 

IHVERSEPERMUTATION OF INTEGER VECTOR 

procedure inversepermb (P ,n) ; value n; integer n; integer array P; 

comment HUCC LIBRARY PROCEDURE NC 06 
Author : J, Boothroyd 

Rearranges the natural number integer elements of P[l:n] 
to effect an inversepermutation 
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This process is best described by :- 

for i :=1 step 1 until n d£ Q[P[i] ] :=i; 
for i:-l step 1 until n do P[i] :=Q[i] 

though NC 06 accomplishes the transformation without the 
use of the auxiliary array Q; 



( 



procedure inversepermb(P,n) ; value n; integer n; integer array P; 
begin ' integer i, j; switch s : -next ; 

for i:=1 step 1 until n do P[i] :=-P[i ] ; 
for n:=n step -1 until 1 do 
begin j:=n; 

next: i :=P[ j] ; if i>0 then begin j:=i; goto next end 
P[j3:=PC-i]; P[-i]:=n 

end 

end i nversepermb ; 



NC 06 
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CRITICAL PATH SCHEDULING 

procedure CRITICALPATH(n,I, J,DIJ,ES,LS,EF,LF 9 TF,FF,i 9 out) ; 
integer i,n; integer array I, J,DIJ,ES ,LS,EF,LF,TF,FF; label out; 
comment HUCC LIBRARY PROCEDURE OPOls 
AUTHOR ACM 74 (improved) : 

Given the total number of jobs n of a project, the pair 
I[k],J[k] representing the kth job as a directed line 
joining event I[k] to event J[k] (I[k]<J[k] ^=1,2,. . . ,n) 
and a duration vector DIJ, the procedure determines the 
earliest starting time ES[k], latest starting time LS[k], 
earliest completion time EF[k], latest completion time LF[k] 
total float TF[k] and free float FF[k], I[l] must be 1 and 
the I[k],J[k] must be in ascending order with no values in 
the sequence 1 to n-1 omitted from the I[k] sequence. The 
critical path is the chain of all jobs with zero total float. 
The following tests are included % 

(a) That I[k]<J[k], (b) All I[k] are in ascending sequence 

(c) No I[k] is missing. 
Failure to meet these requirements causes exit to the 
label parameter out. This exit is so arranged that by 
the use of an actual parameter OUT[i], for example, the 
respective failure paths are (a) 0UT[1], (b) 0UT[2], 
(c) OUT [3]; 



I 
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procedure CRITICALPATH(n, I # J,DIJ,ES,LS f EF,LF,TF f FF,i,out) ; 

integer i,n; integer array I, J f DIJ,ES,LS f EF,LF f TF, FF; label out; 
begin integer k, index, m,tiik,te jk,di jk; integer array ti,te[i:n] ; 
index :=1 ; 

for k:=l step 1 until n do 

begin if ltk]>J[k] then begin i:=1 ; goto out end ; 

if I[k]<index then begin i:=2; goto out end ; 

if I[k]>index and I[k)± index+1 then begin i:=3; goto out end ; 
if I[k]=index+1 then index :=I[k] 

end ; 

for k:=1 step 1 until n do begin ti[k]:=Q; te[k]:=9999 end ; 
for k:=l step 1 until n do 

begin m:=ti[l[k]]+DIJ[k] ; 

if tiEJCk]]<m then ti[j[k)3:=m 
end ti; 
te[J[n]]:=ti[j[n]]; 
for k!=n step -1 until 1 do 

begin m:=te[J[k]]-DIJ[k] ; 

"if te[l[k)]>ro then te[l[k]]:=m 
end te; 
for k:=1 step 1 until n do 

begin tiik:=ti[l£k]]; tejk:=te[J[k]]; di jk:=DIJ[k] ; 
ES[k] : =tiik ; EF[k] : =t iik+di jk ; 
LS[k]:=tejk-dijk; LFfk] :=tejk; 

TT[k] : =te jk-tiik-di jk; FF[k] : =ti [ J[k]]-tiik-di jk 

end 

end CRITIC ALPATH; 



SS01 



REAL RANDOM NUMBER GENERATOR 

real procedure random; 
Comment HUCC LIBRARY PROCEDURE SSOl: 
AUTHOR J. BOOTHROYD : 

A procedure to generate a stream of uniformly distributed 
real random numbers in the interval 4 random 4 1* Random 
integers xr are formed by the operation. 

xr • = ((2^0 - 3) xr - 1) modulo 2 38 

and each real value random is formed from a corresponding 
value of xr by discarding the least significant nine bits. The 
resulting 30 bit truncated integer is then treated as a fraction f 
in the interval < f < 1 and an appropriate adjustment made in 
the integer to real conversion. The period of the sequence is 
thus 2 38 though, by virtue of the truncation, each value of 
random will repeat 512 times in a complete cycle. 

To obtain a uniform distribution in the interval [a,b] use 
the expression a + (b-a) * random. To obtain a symmetrical 
triangular distribution with mean =1 in the range < y < 2 use 
the statement y := random + random. 

Note that random is a procedure without parameters. The 
initial value of random is thus indeterminate. In circumstances 
where this might be inconvenient, as for instance in successive 
program testing runs, it would be preferable to use LIBRARY SS02; 



SS02 



Real Random Number Generator 
real procedure randG2 (x, start); value x, start; Integer x, start 
comment HUCC LIBRARY PROCEDURE SS02: 

AUTHOR J. BOOTHROYD : 

A procedure having the same function as LG1 but with 
facilities for determining the starting value of the 
stream of random numbers. A call of the procedure with 
start = and x » X will cause the first number to be 
generated from 

((2 20 -3) X - 1) modulo 2 38 
Subsequent calls of the form randGZ (0, 1) will generate 
numbers from the value of randG2 at the immediately 
previous call; 



comment LG1 ; 

real procedure random; 

begin own integer xr ; integer three8,co; real f ; switch S:=out ; 

oc:=274877906432; three8:=38; 

elliott(0 ,2, xr f ,0,1,0); 

elliott<5,0,1,0,G,5,xr); 

elliott(5,0, 19,0, 0,4, xr) ; 

elliott(5, 0,18,0,5,7,0); 

elliott(2,0,xr, 0,0,3 ,cc) ; 

elliott<6, 5,4096,0,0, 5, three8); 

elliott(2,0,f,0,4,3,out) ; 
out: random :=f 

end random; 



^Comment LG2; real procedure rand G2(x, start) ; value x, start; integer x, start; 

begin own integer xr,cc, three8; real f ; switch S:=out; 

if start =0 then begin xr:~x; three8:=33; ec: =27487790643 2 end 
elliott (Q rt 2,xr ,0 ,0 , 1 ,0 ) ; 
elliott(5,0,1,G,0,5,xr); 
elliott<5 # 0, 19,0, 0,4, xr); 
elliott<5,0,18,0,5,7,0); 
elliott<2,j,xr,C,0,3,cc); 
elliott(6,5,4096,0,0,5,three8); 
elliott(2,0,f, 0,4,3, out); 
out: rand G2:=rf 

end randG2; 



SS03 



RANDOM NUMBER GENERATOR - PQISSON DISTRIBUTION 

integer procedure poisson(lambda) ; value lambda; real lambda; 
comment HUCC LIBRARY PROCEDURE SS03; 

AUTHOR J. BOOTHROYD i 

A procedure which generates integers x having the 
Poisson probability distribution 

x! 

The numbers are generated from a distribution uniform 
in £ y < 1 by forming x terms of a continued product 
satisfying 

71 Y2 Y3 Y x < e X 

For large values of X the method is therefore not 
particularly efficient, and it may be preferable to use 
the fact that for large X the Poisson distribution 
tends to a normal distribution and employ SS04; 



integer procedure poisson(lambda) ; value lambda; real lambda; 
begin ovm integer xr; integer threes, n; real elambda,try; switch 

threes : = 38; n:s-1; try: = 1.Q; elambda:=exp( -lambda) ; 
L: elliott(o,2,xr, 0,0, 1,3) J 

elliott(5,0, 1 ,0, 3, 5,xr) ; 

elliott<5,0, 19, l,0,4,xr) ; 

elliott(5,0, 18, ,5,7,3); 

elliott(l,6,xr,0,5,4,29>; 

elliott(5, 5,9,0,6,5,4096); 

elliott(0,5,three8, 0,6,3, try); 

elliott(2,0,try,0 f 6,2, elambda) ; 

elliott(2,2,n f 3,4,1, L) ; 

elliott<7,3,4,1,4,3,1); 
poisson:= n 

endp oisson; 
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RANDOM NUMBER GENERATOR - STANDARD NORMAL DISTRIBUTION . 

real procedure normal ; 
comment HUCC LIBRARY PROCEDURE SS04: 
AUTHOR J. BOOTHROYD : 

Generates random numbers with a standard normal distribution 

(mean zero and unit standard deviation). The numbers are 

generated from a distribution uniform in < y < 1 

by summing twelve values of y and adjusting the mean 

by subtracting 6,0 

12 

normal := £y. - 6.0 
1 1 

To generate numbers having a mean u and standard deviation 
a use the transformation 

musigma := sigma*normal + mu ; 



1 

j SS04 

real procedure normal; 

begin own integer xr ; integer three8, count ; real sum; switch s:= L; 
threes : =38; count := -11; sum:= -6.0; 
L: elliott(o,2,xr,0,0,1,0); 
elliott<5,0,1,0,0,5,xr) ; 
elliott(5, 0,19,0,0, 4, xr) ; 
elliott(5, 0, 18,0,5,7,0) ; 
elliott(l,6,xr,0,5,4,29) ; 
elliott<5,5,9,0,6,5,4096); 
elliott( 3,5, three8 ,0,6,0, sum) ; 
elliott(2,o » sum, 0, 3, 2, count) ; 
elliott(4,1,L,0,0,0,0); 
elliott(7,3,4,1,4,3,1); 
normal := sum 

end normal ; 



SS05 



GENERATE LARGE INTEGER RANDOM NUMBERS 



integer procedure bigrn; 
comment HUCC LIBRARY PROCEDURE SS05: 
AUTHOR J* BOOTHROYD s 



A procedure to generate, equiprobably , integers in 
the range 4 bigrn < 2 38 -l. There is thus a 0.5 
probability that bigrn > 2 37 . The numbers are formed 
by the recurrence 

xr: = ((2 20 -3)xr~l)modulo 2 38 . 
For repeatable program testing purposes it may be 
preferable to use SS06; 

SS06 



GENERATE LARGE INTEGER RANDOM NUMBERS 



integer procedure bigrnx(x, start) ; value x, start; integer x, start; 
comment HUCC LIBRARY PROCEDURE SS06: 
AUTHOR J. BOOTHROYD : 

Performs the same function as SS05 but with facilities 
for initialising the starting value of the stream of random 
numbers. With startO the value of bigrnx is generated 
from the value supplied by the parameter x. For start=l 
the value of x is ignored and operation of SS06 is the same 
as SS05; 



SS07 



GENERATE SMALL INTEGER RANDOM NUMBERS 

integer procedure rn (n) ; value n ; integer n ; 
comment HUCC LIBRARY PROCEDURE SS07: 
AUTHOR J, BOOTHROYD : 

A procedure to generate, equiprobably , integers in 
the range 4 m < 2 n . The basic generation scheme 
is the recurrence 

xr: = ((2 20 -3)xr-l)modulo 2 38 
The value given by the top n bits of each generated 
xr is assigned to m. Where repeatable program testing 
results are required, it may be preferable to use SS08; 



SS08 



GENERATE SMALL INTEGER RANDOM NUMBERS 

integer procedure m(x,s tart,n) ; value x, start, n; integer x, start, n; 
comment HUCC LIBRARY PROCEDURE SS08: 
AUTHOR J. BOOTHROYD : 

A procedure performing the same function as SS07 
but with facilities for determining the starting 
value of the stream of random numbers. With 
start=0 the value of the first generated number is 
determined by the value of x. For start=l, x is 
ignored, and the procedure is identical with SS07 ; 



integer procedure bigrn; 
begin own integer xr; 

elliott(0,2,xr,0,0,1,0); 

elliott(5,C, 1,0,0, 5,xr) ; 

elliott(5,0, 19,0, 0,4, xr) j 

elliott(5,Q, 18,0,5,7,0) ; 

elliott< 2, 0,xr, 0,^,0,0); 

elliott<7,3,4,1,4,3,1); 

bigrn := xr - } 

endbigrn; i ss06 



integer procedure bigrnx(x, start) ; value x, start ; integer x, start; 
begin own integer xr; 

if start=0 then xr: = x; 

elliott(0,2,xr,0,0, 1,0) ; 

elliot t (5,0,1,0,0,5, xr ) j 

elliott(5,0,19,0,0,4,xr) j 

elliott(5,0,18,0,5,7,0) ; 

elliott(2,0,xr, 0,0,0,0) ; 

elliott(7, 3,4,1,4,3, 1); 
bigrnx:= xr 

end bigrnx; 
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integer procedure rn(n) ; value n; integer n; 
begin own integer xr; 

" elliott(0, 2, xr, 0,0,1,0); 

elliott(5, 0,1,0,0, 5,xr) ; 

elliott(5,0, 19, 0,0, 4, xr) ; 

elliott( 5,0, 18,0,5,7,0); 

elliott(1,6,xr,0,6,7,n) ; 

elliott(5, 4,0,0,2, 0,n) ; 

elliott (7, 3,4, 1,4,3,1); 
rn:=n 

endrn; 
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integer procedure rnx(x, start, n) ; value x, start, n; integer x, start, n; 
begin own integer xr; 

if start=0 then xr: = x; 

elliott (C, 2, xr,0, 0,1,0) ; 

elliott(5,0, 1,0,n,5,xr) ; 

elliott(5,0,19,0,0,4,xr); 

elliott(5, 0, 18,0 , 5, 7,0) ; 

elliott (1, 6, xr, 0,6, 7, n) ; 

elliott(5,4,0,0,2,0,n) ; 

elliott(7,3,4,1,4,3, 1) ; 
rnx: =n 

endrnx; 



ZM01 

OUTPUT BINARY VALUE 

procedure binary (i s g) ; 
value i,g; integer i,g; 
comment HUCC LIBRARY PROCEDURE ZMOlr 
AUTHOR: W. G« WARNEt ' 

Prints on current device the binary value of the 
integer i in groups of g digits separated by a 
space starting from the most significant binary digit. 
No character is output before the first digit. 
If g « no spaces are output to cause grouping. 



ZM02 

OUTPUT OCTAL VALUE 

procedure octal(i f g) j 
value i,g; integer i,g; 
comment HUCC LIBRARY PROCEDURE ZM02 : 
AUTHOR: W. G. WARNE: 

Outputs on current device the octal value of the 
integer i in groups of g octal digits, separated 
by a space, starting from the most significant octal 
digit. No character is output before the first 
octal digit* 

If g * the integer is output as a single group. 



< 



p rocedure binary< i)in groups of :(g) ; 
valte ifgr; integ er i,g; 
begin integer n,gc; 

switch sw := neg, cont; 

gc:=0 ; 

lor n:=1 steg 1 until 39 do 

begin elliott(3 ,0 ,i ,0 ,4,1 ,neg) ; 
print £0?; 
goto cont; 

neg: £rint £1?; 

cont; gc:=gc+1 ; 

if gc=g then begin print ££s?? 

elliot t<3 ,0 , i ,0 , 5 , 5 ,1 ) ; 

el liott< 7 , 3 ,4 ,1 . ,4 , 3 ,1 ) ; 

elliott(2,0,i,0,0 ,0 ,0) ; 

end 



er^\ binary; 



gc:=0 end; 



procedure octal(i) in groups of :<g) ; 

valtia integ e r i a g; 

begin int eger n,gc , temp ,seven; 

gcs=:0 ; 

seven: =7; 

for n:=36 step -3 until do 

begin elliott<3 J* t i s 0»6 t 7,n) ; 

elliott<f>,1 ,0 & J) B 3 ,eeveii) ; 

elliott( 2 ,0 ^temp/I »2 »2 ^gc) ; 

print sameline ,special(1 J # digi ts(1 > ,temp ; 

if gc=g the n b egin p r i n t ££s??; gci=0 end 

end 

ena rctal; 



ZM02 
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SENSE NUMBER GENERATOR KEY 

boolean procedure ng(i) ; value i; Integer 1; 
comment HUCC LIBRARY PROCEDURE ZM03s 
AUTHOR ELLIOTT AUTOMATION: 

This procedure, when called as ng(x) where x is 
an integer, 1 < x < 39, takes the value true if the 
key numbered x on the word generator is depressed 
at the time of call. If key x is not depressed, the 
value of ng is false. A typical call would then be 
if ng(x) then ; 



ZM04 



INPUT NEWLINE STRINGS 

procedure charead( array) ; integer array array; 
comment HUCC LIBRARY PROCEDURE ZM04: 
AUTHOR P. W. FORD : 

A call of charead causes the input tape to be searched 
for the first occurrence of a non newline non blank 
character immediately following a newline 
i.e. blank tape and a succession of newlines is ignored. 
This first character and all subsequent non blank 
characters up to but not including the next newline are 
stored in the integer array array in a form consistent 
with subsequent use of the procedure outs tring. The 
character string is stored starting in array [1] and 
adequate bounds for this array must be declared (see 503 
TECH MANUAL 2.1.3,2 section 2.6). There is no protection 
on the bounds of array and the input of a string longer 
than can be accommodated is likely to be discovered by a 
subscript overflow failure when outs tring is called; 
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boolean pr ocedu ro ng( i) ; value i ; integer i • 
begin sv/itch s:=l; 

elliott(0 ,2,0,0,0 ,0 ,0) ; 
elliott(0 ,0 , i ,1 ,5 ,4 f 81 91 ) ; 
elliottd ,6,i,0,7,0 ,0) ; 
elliott( 2 , 3 , i ,0 ,4 , 3 , l) ; 
1; ng:=i $ 
end ng(i) ; 
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procedure chare ad (array) ; 
integer array array: 

begin integer i f j,k,l,m; 

switch ss:=search,next, read, store, shift # test f out , back, on; 
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end charead; 



ZM05 
ZM06 

Algol A«D,T, Read and Punch Binary 

Two procedures designed to meet the requirements of programs 
using paper tape as temporary storage so that partial results output 
by one program are used solely as subsequent input data for the same 
or another program. In these cases the most appropriate information 
format is pure binary and further efficiency has been obtained by 
taking advantage of the interrupt facilities of the 503 so that input 
and output takes place simultaneously with computation or the use of 
peripheral devices other than those used by ZM05 or ZM06, which may 
operate together. 

Both procedures use a boolean parameter flag to indicate when 
a transfer initiated by a call of ZM05 or ZM06 is complete. Programs 
using these procedures should not access the data transfer locations 
unless flag = true . 

Realignment of a tape, output by ZM05, at subsequent use by 
ZM06 may be accomplished in several ways:- 

(a) By arranging to output a correctly terminated number 
which is read in and checked before a call of ZM06 

(b) By arranging to output an identifiable character and, 
on reinput, to execute a buffer search for the same 
character, or more simply but without any check, execute 
a call of the procedure advance which ignores blank tape 
and erase characters, 

ZM05 

procedure out(a,n,f lag) ; value n; integer n; array a; boolean flag; 

comment HUCC LIBRARY PROCEDURE ZM05 
AUTHOR : W.G. WARNE 

A procedure which outputs, in binary* on punch 1, the first 
n(>l) locations of array a, i.e. by rows. The parameter flag 
is assigned the value true when the transfer is complete. 
Total time is approximately 60n milliseconds of which 99% 
is available for simultaneous computation; 



ZM06 

procedure in(a,n,f lag) ; value n; integer n; array a; boolean flag; 

comment HUCC LIBRARY PROCEDURE ZM06 
AUTHOR : W.G. WARNE 

A procedure which reads, from reader 1, tapes previously 
output by ZM05. (see the notes above corocsrning 
realignment of tape). Data from the tape is input to 
the first n(:>l) locations of array a. Total time is 
approximately 6n milliseconds of which 90% is available 
for simultaneous computation; 



ZM05 

procedure out(a,n,f lag) ; 
value n; 
array a; 
boolean flag; 
integer n; 

begin own integer aa, i,ab, temp, hold, x, mask, imask; 

switch s: = f in, LI ,L2,L5; 
mask := 127; 
x := 32; 
imask := 223; 
flag := false ; 
aa := address(a); 
ab : = aa + n ~ 1 ; 
goto L1 ; 
1*2: i := i - 262144; 

elliott (3, 0,7886,0, 0,3, imask); 

elliott(0, 4, x, 0,2, 0,7886); 

elliott(6, 7, 7886,0, 7,2,0); 

elliott(0,2,i,0,5,5,20); 

elliott (2, 6, 7895,0, 1 ,0,7892) ; 

elliott (2 ,0 ,hold ,0,2,6, 7893) ; 

elliott(2, 6, 7894, 0,3,0, fin); 

elliott(1,6,4,0,0,2,0); 

elliott (5 ,5 , 3,0,2,0, 7896) ; 

for aa := aa step 1 until ab do 

begin elliott(0,6,aa,1 ,0,4,0); 

elliott (2,0, temp, 0,0, 0,0); 

for i := 35 step -7 until do 

begin elliott(3,0,ten$>,0,6,7,i); 

elliott (5, 1,0,0,0,3, mask); 

elliott(2,0,x,1 ,7,4,0); 

elliott(6, 6, 7893, 0,0,0,0); 
LI: elliott<7,3,i,0,4,0,L2); 

elliott(7,5,7168,0,5,5,2); 

elliott(7,3,x, 1,4,3,1); 

elliott(4,1,L5,0,0,C,0); 

end ; 

end ; 

flag := true ; 

elliott(3, 0, hold, 0,2, 0,7892) ; 

elliott(3,0, imask, 0,0, 3, 7886); 

elliott (2,0, 7886 ,1,7,2,0); 
L5: elliott(6, 6, 7893,0,0, 0,0); ^ 

fin: end ; * 



procedure in (a, n, flag); 
value n; integer n; 
array a; 
boolean flag; 

begin own integer aa,i,ab, temp, hold, five, imask; 
switch s:= f in, LI ,L2,L3,L4,L5; 
flag:= false ; 
aa:=address(a); 
ab:=aa+n-1 ; 
five := 5; 
temp := 64; 
imask := 191; 
goto LI; 
L2: i:=i-262144; 

elliott(3, 0,7886, 0,0, 3, imask); 
elliott( 0,4, temp, 0,2, 0,7886); 
elliott(6, 7,7886, 0,7, 2,0) ; 
elliott(0,2,i,0,5,5,20); 
elliott<2, 6, 7888, 0,1, 0,7887); 
ell iot t ( 2 , C , hold ,0,2,6, 7889 ) ; 
elliott(3, 0, five, 0,2,1 ,i) ; 
elliott(2,6,7890,0,3,0,fin); 
elliott(1,6,4,0,0,2,0); 
elliott(5,5,3,G,2,0,7891); 
elliott(0,6,0,0,4,0,L3); 
L4: elliott(2, 6, temp, 0,3,0, five); 

elliott(2,1,i,0,6,6,7888); 
M: elliott(7,3,i,0,4,0,L2); 

elliott(7,5,7168,0,4,1,L5); 
elliott(3,0,temp,G,5,5,7); 
L3: elliott(7, 1 ,0,o ,2,0, temp) ; 

elliott<3,2,i,0,4,1,L5); 
elliott(3,0,temp,0,6,7,aa) ; 
elliott(2,0,0,0,3,2,aa); 
elliott(0,5,ab,0,4, 1 ,L4) ; 
flag:=true; 

elliott(3,0, hold, 0,2, 0,7887); 

elliott(3,G, imask, 0,0,3,7886); 

elliott(2, 0,7886, 1,7, 2,0); 
L5; elliott<6,6,7888,0,0,0,0); 
fin: end; 



ZM07 

COPY LEGEND FROM DATA TAPE 

procedure copy(n) ; value n; integer n; 
comment HUCC LIBRARY PROCEDURE ZM07: 
AUTHOR : W.6. WABNE: 

Copies characters from reader (1) to punch (1) up to but not 
including the first character having binary value n where 
5l n <. k27 in 8-hole mode or £ n <. 31 in 5-hole mode; 



ZM08 

Read ALGOL BUFFER 

integer procedure BUFFER (d); value d; integer d; 

comment KUCC LIBRARY PROCEDURE ZM08: 
AUTHOR : W.G. HARNE: 

Takes the value of the character in buffer(d) where d(l,2 or 
is the device number. If d is not 1,2 or 3 the message 
"BUFFER error 11 is displayed followed by a data wait. 
Continuation assigns the value zero to BUFFER, The 
procedure has been tested with Algol 1/3 Tapes 1 and 2 only 
on a basic machine. There is thus no guarantee that it is 
operative with the backing store version of Algol 1; 



ZM07 



procedure copy(n) ; 
value n; integer n; 

comment Copies characters from reader 1 to punch 1 until a character with 

binary value n is met. This character is not copied. 

Note 0<n<l27 in 8 hole mode and <Kn<31 in 5 hole mode.; 
begin integer t; 

switch L := L1 f L2; 
LI : elliott<0,6,0, 0,7,1,0) ; 

elliott(2,0,t,0,0,5,n) ; 

elliott(4,6,L2,0,6,7,t) ; 

elliott<7,4,G,0, 4,4,L1> ; 

L2: 
end copy; 



ZM08 



integer procedure BUFFER(d) ; 

comment Takes binary value of character in buffer for device d (d=1 , 2 or 3).; 
value d; integer d; 
if d < 1 or d > 3 then 

begin print punch(3), ££1?BUFFER error?, wait; 
BUFFER := 

end 

else begin elliott(3, 0,7920,0,2, 4,d) ; 

elliott(0,6,d,1 ,0,4,79) ; 
elliott<2,0,d,0,0 f 0,0) ; 
BUFFER := d 
end BUFFER; 
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PRINT MATRIX 

procedure matprint(fields , rows , cols ,i ,j ,margin ,gap , var) ; 

value fields, rows, cols; integer f ields , rows , cols ,i ,j ; real var; 

string margin, gap; 

comment HUCC LIBRARY PROCEDURE ZP02 : 
AUTHORS J. BOOTHROYD S 



A Jensen procedure for printing the elements of either one 
or two dimensional arrays with facilities for specifying 

(a) fields .... the number of numbers printed per line 

(b) margin .... the prefix required at the start of each line 

(c) gap the number of spaces (or other symbols) 

printed between numbers. 
Notes 1. i and j are the Jensen parameters, with ranges 
i, 1(1) rows, j ,1(1) cols. 

2. var must be called as a real variable, A[i,j] 

for a matrix A[lsrcws,l:cols] or a[j] for a vector 
a[lscols] . For arrays with other than unit lower 
bound the appropriate adjustments should be made 
to the subscript expressions and to the values used 
for rows and cols. For example, a in the case 
of TABLE [0:6,0: 10] the appropriate parameters 
would be 

rows = 7 cols = 11 

var - TABLE [i-l/j-1] 

3. The procedure does not specify the style of number 
printing. The user has the option of any of the 
following methods of determining the required format; 

(a) By calling the required setting procedure 
( aligned (m,n) , digits (n)etc) globally. 

(b) By calling the required setting procedure 
locally in a print statement. 

e.g. print aligned(2,3) ,matprint( ) ; 



_2- ZPG2 cont. 

3 (cont.) 

(c) By declaring a private format procedure to 

permit optional setting procedures depending 
on circumstances 

e.g. real procedure format(x) ; value x; real x; 

begin if abs(x)<0.1 then scaled (7) else aligned(5,6) 
format t~ x 

end ; 

Within matprint this procedure would be used as 
an actual parame ter as in 

matprint(5 ,10 ,10,i , j ,££s5?? ,££s4?? f format(A[i f j]» 
4. It is the users responsibility to ensure that the total 

character requirements of the setting procedure, margin and 
gap are consistent with the requested value of fields and the 
limit of 120 characters per line. 

Examples ; Tc print A[l:20,l:20] with 5 spaces as margin, 

3 spaces between numbers and 6 numbers per line. 
matprint(6,20 9 20,i,j,££s5??,££s3??,A[i,j]) 

To print LIST[l:n] with each line starting one 
tab from the paper margin, 7 numbers per line, 

4 spaces between numbers • 
matprint(7,l,n,i,j,££t??,££s4??,LIST[j]) 



r 



ZP02 



procedure matpr int (fields , rows , cols , i , j , margin, gap , var) ; value fields, rows, cols; 
integer fields, rows, cols, i,j ; real var; string margin, gap; 
begin integer k,line,np; switch s : =mor e ; 
for i:=l step 1 until rows do 

begin print ££12??; j:=0; np:=cols; 
more: np:=np-fields; line;=f ields; 

if np < then line :=np+f ields; 
print margin; sameline; 
for k:=l step 1 until line do 
begin j:=j+l ; 

if k i 1 then print gap; 
print var 

end ; 

print ££1??; if np > 3 then goto more 

end i 

end matjrint ; 
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CONTROL FLEXOWRITER PAGE PRINTING 

procedure pagequt(string f type,resultline f array, start) ; 
value type, resultline, start; string string; 
integer type, resultline, start; integer array array; 
comment HUCC LIBRARY PROCEDURE ZP03: 
AUTHOR P. W. FORD 

A procedure to control the automatic pageing of results 
printed on Lams on Paragon Paraflo paper (Form No. 1112). 
The number of resultlines is internally set at 50 and 
change to a new page occurs if the next batch of output 
would cause this limit to be exceeded. Page numbering 
is provided and facilities are provided for including on 
each page indicative information such as job identification, 
result headings, etc. A call of page out must occur before 
each print statement or group of print statements; 

Parameters 

t yP e an integer with permissible values 0,1,2. 

type * prepares the procedure by initialising the 

the line and page counters, 
type = 1 for normal page control. 

type = 2 initiates an immediate page change and resets 
the page counter. 

resultline an integer specifying the number of lines of output in 
the immediately following print statement (s) . This 
parameter is ignored, and may be 0, if type = or 2. 
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Parameters (cont.) 

array the identifier of an array containing some job or 

other indicative string to be output at the head of 
each newpage on the line above the page number. If 
this facility is not required (start « 0, see below) 
use any declared integer array identifier if one exists 
otherwise declare a minimal array, e.g. A[l :1] for 
this purpose. 

start an integer parameter which controls the outstring 

facility used for job identification. If start = 
the outstring facility is ignored, while start £ 
causes the string starting at array [start] to be output 
at the head of each newpage. Note that start may be an 
integer constant, and does not have to be an integer 
variable identifier. 

string Array string (including the empty string £?) . This 

is output on each newpage following the line bearing 
the page number. This facility permits the printing 
of result headings , etc, 

NOTES STRINGS should NOT contain newline characters. 

This procedure does not control line changing. This 
must be done in the print statements whi ch pageout 
monitors. The procedure must be initialised once 
with type « 0, and this is preferably done in the 
housekeeping. 



procedure pageout( string, type, resultline,array # start) ; 
value type, result line, start ; 
string string; 

integer type, result line, start; 
integer array array; 
begin own integer line, page; 
switch s8s:=LT 9 L2,L3| 
if type=2 then goto LI ; 
if type=:1 then 

begin if line+resultline>50 then 

LI: begin for line:=line+l while line<62 do print ££!?? 
if type=2 then goto L3; 
goto L2 

end 

else 1 i ne : =1 i ne+r esul t 1 i ne 

end; 

if type=0 then 

L3 : begin page t =1 ; 

L2; print £fil??; 

if start * then outstring( array , start) ; 

grint fi£lt8?page?,sameline,digitsf3),page,££12?? r string7 
pagej=page+1 ; 

liner=if type=t then resultline else f) 
end — — 
end procedure pageout; 



) 



